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ABSTBiiCT 

In the present study mathematical models for the blast 
furnace process have been reviewed. A decoupled multivariate 
stochastic model has been developed using multivariate time 
series analysis. This method is a generalization of the method 
proposed earlier by Phadke et al. 

Data have been collected for the blast furnace at Bokaro -J 
Steel Limited, Bokaro Steel City, The input varialsles considered 
in the analysis were sinter-to-coke ratio, blast humidity and 
blast flow rate. The output variables taken into account were 
hot metal temperature, silicon and sulphur contents of hot 
metal. The original series were first 'prewhitened’ by means of 
unvariate modelling. The prewhitened series which are normally 
distributed and serially independent were used in multivariate 
time series analysis. 

The method uses Gram-Schmidt procedure to obtain orthogonal 
vectors. A principal component model was derived in terms of 
prewhitened series and orthogonal vectors. A multivariate model 
was then developed using cross-correlation function between 
orthogonal vectors and multiple regression analysis. The 
orthogonal vectors were found to be serially and mutually un- 
correlated andhence the multivariate model was also a principal 
component model. 
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Using this multivariate model, transfer function models 
were developed, both in terms of the prewhitened input-output 
variables; and in terms of actual input-output variables by 
recoupling the multivariate model. It was found that the 
input-output relationship consists of terms involving shift 
operator which represent the time delay and other terms which 
represent the time constant. The average time lag between each 
input e.nd output variables was found to be of the order of 
0.85 to 1.10 cast intervals. 

These results are generally in agreement with earlier 
results. But they indicate greater details of the transfor- 
mation relationships. Particularly important are the direct 
interactions of the order of 2 to 3 cast intervals between the 
input and output variables. However no feedback control of 
the order of one cast interval or more were indicated in the 
study. 
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CHAPTER 1 


INTRODTJCTIOH 

1.1 aSHSMXj 

Blast furnace is a major equipment of an integrated 
iron and steel plant which produces pig iron or hot metal of 
desired quality. The principal raw materials used in the 
blast furnace are iron ore, sinter, coke, limestone and manga- 
nese ore. Other materials such as quartzite, dolomite, scrap, 
open hearth slag etc are used depending upon the requirement 
and their availability. 

Fundamentally, the blast furnace is a counter-current 
apparatus in which chemical and heat transfers take place. 

The solids added at the top of the furnace are usually coke 
and sinter. Other materials are added to take care of fluc- 
tuations in the sinter composition and/or to get the desired 
composition of the metal. Sinter is an agglomerate of ore 
fines, limestone and some essential additives which facilitate 
the production of iron. The gaseous phase is a mixture of 
CO, CO 2 , H 2 , H 2 O and 1^2 which is more and more oxidized as it 
ascends the furnace and releases heat. The initial gas 
constituents are CO, H 2 and IT 2 high temperature. They 
result from combustion of a part of coke by the blast air and 
decomposition of steam injected along with the blast. The 
blast air is preheated to about 1000°C. The gas heats the 
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solid burden, reduces iron oxides to iron which, is finally 
melted and collected in the hearth. The slag which is formed 
with the gangue material floats on molten pig iron. These 
liquid products are tapped after every two to three hours. 

The blast furnace is characterized by a considerable 
number of inputs and outputs which together define the operating 
status of the furnace. The inputs or independent variables 
are defined as those which affect the behaviour of the system. 
These variables include those which can be directly manipulated 
by the operator (manipula table variables or controllable 
variables) as well as those which are beyond his direct control 
(disturbances or noise in the system). Outputs or dependent 
variables are the remaining variables of the process which are 
the result of the input variables and are needed to describe 
th^ operating conditions of the process. In other words these 
are responses of the system to the inputs. The manipulatable 
variables may be classified into two principal categories: 
charge variables such as sinter-to-coke ratio; and operating 
variables such as blast temperature, blast humidity and blast 
flow rate. The smooth operation of the blast furnace is affected 
by many disturbing factors. The significant disturbances are 
chemical and physical composition of the charge material, 
porosity of the bed etc. The outputs or dependent variables 
are the temperature of hot metal and its sulphur and silicon 
content. 
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1.^ EEBD FOR A MODEL; 

The advent of fast and sophisticated I.D. steelmaking 
process demands great improvements in the control of iron- 
making process and in particular better control of silicon, and 
sulphur content of the molten metal. Models of the furnace 
dynamics are needed to design a control scheme to minimize 
variations in the silicon content in pig iron from cast to 
cast, A well regulated silicon content is likely to lead to 
smooth and economical ironmaking as well as steelmaking proc esses c 
The ©Gihststentnsilicon content implies stable thermal conditions 
in the lower part of the furnace. Studies of control systems 
are based on a mathematical model of the process involved, A 
large number of blast furnace models have been proposed in 
past ranging from completely theoretical to fully empirical 
ones, A review has been published earlier [l]. 

Blast furnace contains a number of inputs and outputs 
and in order to achieve optimal control, it is necessary to 
identify (i) the transfer function relating each manipulatable 
variable with each output and (ii) the disturbance or noise 
present in the system. The classical methods of estimating 
the transfer functions are based on studying the result of 
specific deterministic perturbations of an input while maintain- 
ing other inputs at their steady levels. The perturbation may 
be of the form of a step, pulse or sinusoid. However, such 
procedures have not a,lways been successful. This is because 
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for perturbations of a magnitude that are relevant and tolera- 
ble, the response of the system may be masked by uncontrollable 
disturbances in the system, and the large capacity of the 
system may absorb permissible magnitudes of the perturbation. 
Further other inputs have to be maintained at their steady 
levels which is very difficult in the plant environment. 

The sequence of values of a variable observed along 
time is called a time series. The sequencesof values of one 
variable constitutes a univariate time series. The sequences 
of values of number of variables constitute a multivariate time 
series. The statistical analysis of multivariate time series 
takes into account the disturbances present in the system. 
Multivariate modelling involves the representation of the 
relationships among input and output variables and is consi- 
dered in this study. 

1.3 OBJECTIVES OF THE STUDY: 

The objectives of the present study are as follows; 

(a) to fit a s-uitable univariate model to each of the 
variables under consideration so as to obtain residual series 
which are pure random. If the residuals have a normal distri- 
bution, they may be called as ’prewhitened series*. 

(b) to fit a suitable multivariate model to the serially 
independent residual series in terms of serially and mutually 
independent random components, and 
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(c) to validate the fitted multivariate model and 
hence derive hy suitable matrix operations, a transfer function 
model from the multivariate model. 

This multivariate model is called a decoupled model 
because, all series are first prewhitened by univariate time 
series analysis. This approach for modelling multiple time 
series is advantageous because; 

(i) as the parameters of univariate and multivariate 
models are estimated separately, their numbers are comparable 
to those of univariate and multivariate models and fhe diffi- 
culties involved in simultaneously estimating all parameters 
are avoided; and 

(ii) the most appropriate form and order of model can be ’ 
chosen independently for each variable, and hence there is 
greater flexibility in the choice of the models. 

1.4 ORGAHIZITION OP THE STUDY; 

The study is reported in the following sequence; 

(i) Chapter 2 is a review of mathematical models that 
have been proposed in the past, 

(ii) The technique of univariate time series analysis 
are described in Chapter 3' in the following sequence - Eirst 
the general time series models are described and then the 
similarity between the linear difference equation and linear 
differential equation is given. The steps for fitting univariate 
time series models are then described. These include 
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identification of a suitable model, estimation of the para- 
meters of the identified model and the diagnostic checking of 
the residua,ls for their normality and serial independence. 

They are applied to field data from Bokaro Steel Plant. 

(iii) In Chapter 4, a multivariate model is fitted 
to the serially uncorrelated residuals of the unvariate models. 
The method proposed by Phadke et al. [2] has been generalized 
and used in the modelling of the process. The transfer fxmction 
model is then developed from the decoupled multivariate time 
series model. 

(iv) Finally, in Chapter 5 conclusions are drawn on 
the basis of results of the present study and some suggestions 
for future work are given. 
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CHAPTER 2 


MATHEMiTICAL MODELS OP BLAST PURHACE PROCESS 


The mathematical models for blast furnace process may 
broadly be classified as deterministic or probabilistic models. 

A general classification of the models is given in Table 2.1. 

2.1 DETERMINISTIC MODELS: 

A deterministic model is one where the future behaviour 
of a physical system is expressed uniquely by a set of algebraic, 
differential or integro-differential equations and hence there 
is no uncertainty or randomness in the future outcomes. Deter- 
ministic models may further be subdivided into two groups- steady 
state and dynamic models. 

2.1.1. Steady- St ate Models : 

Steady state models deal with equilibrium conditions of 
a process. Considerable amount of work has been done in past 
on formulation of mathematical models for blast furnace assuming 
steady state conditions, i.e. assuming that the operating 
conditions of the furnace do not change with time. These models 
are based on mass and energy balances and may be divided into 
three groups: 

1. Models based on overall heat- and material balances 

2. Models based on stagewise heat- and material balances 
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3. Models based on differential heat and material 
balances 

M odels Based on Overall Heat and Material Balances ; Overall 
balances are commonly used to check the consistency of the 
plant data. Models belonging to this group include those of 
Joseph et al.[3j, Marshall [4], Dancy et al, [5], Lander et al. 
[6-7], Pazzalari et al [8], Khromov et al [9], lahiri et al.[lO], 
and Pokhvisnev et al [11-13] and have been siommarized in Table 
2 . 2 . 

Joseph et al [3] have derived carbon balance, oxygen 
balance and heat balance from the known operating conditions 
and the average top gas analysis. Marshall [4], related blast 
furnace operating variables with coke rate and production rate 
with the help of carbon balance, heat balance, wind rate and 
coke burnt at the tuyeres. Dancy et al [5] calculated -the 
production rate and coke rate for a blast furnace with oxygen 
enrichment a.nd steam injection through enthalpy and oxygen 
balances. It was assumed that the amount of heat required to 
produce a unit of iron from a given burden is constant irrespective 
of blast composition, provided that the burden composition, slag 
composition and the stone ratio do not alter, thimigh the 
reactivity of the raw material, the burden size, porosity, flow 
pattern all have effect on the amount of energy required to 
reduce the ore to iron. The model was applicable to changes 
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brought about blast modification only. Lander et al [6,7] 

modified Dancy’s approach which enabled them to calcula,te the 
effect of variations in the burden composition on the production 
rate and the coke rate. This was accomplished by defining a 
quantity known as furnace characteristic constant obtained from 
reference operating conditions. The composition and temperature 
of top gas were assumed to be similar to that of reference 
conditions. The model was applicable only to small perturbations 
in the operating conditions from that of the reference conditions. 
Fazzalari et al. [8] developed a computer model which consisted 
of. overall carbon balance, heat bale^nce, absolute top gas 
temperature, theoretical flame temperature , production rate, coke 
rate and shaft efficiency. Khromov et al, [9] developed a 
model for stabilization of blast furnace based on a variable 
defined as heat unbalance , that is caused by previous charging 
cycles. The heat unbalance was calculated as the difference 
between the heat required and the heat available in the furnace 
per kg of the carbon in the dry top gas divided by the heat 
available. Lahiri et al. [10] developed a model for calculating 
coke rate and optimum gasification rate utilizing the principle 
of optimum gasification. He assumed that carbon monoxide is 
mainly supplied by the exothermic reaction C + 4’02 = CO and only 
the deficient amount is supplied by the endothermic gasification 
reaction G + CO 2 = 2G0. The influence of charge composition 
and other parameters on coke rate were studied. The. effects of 
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gSngue in ore, of ash in the coke and of silica in the flux 
on the coke rate as predicted by the model were not as signi- 
ficant as in the actua.1 practice. Pokhvisnev et al.[ll-13] 

CT C 

developed two thermal state indices and M for monitoring 
thermal state of the blast furnace. Indices and re- 
presented the heat input in the blast furnace per unit of the 

oxygen of the charge gasified. was calculated from the 

c ii 

top gas composition and M was calculated from the composition 
of the charge, blast and injected additives. Both indices 
had positive correlation with silica content of pig iron. 

Models Based on Stagewise Heat and Material Balances ; Models 
which belong to this group were presented by Reichardt [14], 
Ridgion [15], Hodge [16-18], Pokhvisnev et al.[l9], Dovgalyuk 
et al. [20], ¥artme.nn [21], Staib et al. [22-24] , Daconsine[25] , 
VanLangen et al. [26-30], Flierman et al. [31] and Rist et al.[32] 
and have been summarised in Table 2.3. 

Reichardt [14] studied the longitudinal distribution of 
temperature of gas and solid particles in the furnace and 
demonstrated the existence of the ’thermal pinch point' which 
is associated with the sudden increase in the overall co.pacity 
of the charge. The thermal pinch point is defined as a point 
where the solids and gaseous phases have the same temperature. 
Ridgion [15] developed a model for stagewise heat balance 
assuming that the temperature of the gas at any level in the 





1 AIM OR QDHTROL OBJECTIVE NO, 

NO.OF :BEGiaNS 

1 ccutrolabie 

_ - r - 2 " ’ ’ ’ ’ 5 ' ' ' 

,:J2, .1 



Stagemse -heat ’balance Heat requirement of spli<i 
(heat ■rhfluirei. for ‘ solid phase is calqula-ted in 5GC>G 
and gaseo^ phas%) , steps from to ISOO^C,' 

Reactions fccurring 
divided in '3 groups;- ® , 
Relating to burden cons- 
.titueniB other than iron 
compnmds^ Iron reduction, 
and sola loss* 


Thcee:- (1) Removal of 
moisture (2) Preheating, 
calcination, indirect 
reduction, direct redu- 
ction j reduction with H 2 
(3) Direct reduction, 
slag formation, liauification 
tion and solution of ele- 
ments, high temp •combustion. 


Coke rate and 
Production rate 


Automatic control of the 
it-ral. thermal , state of blast 
furnace, ^ 


Two;- Direct reduction 

(Lower Zone) , and indirect 

reduction (TJfjper Zone) . cr propcrl 

hatu'rl'. g; 
lorer rcnc 

temp, , bli 
humidity, i 


b.al. Computer control of 
the thermal state of 
the blast furnace 


Two;- Direct reduction Ore to - : 
and indirect reduction Coke rati 
zones. 


AMo : 


M&M 


EQUATIONS I DMSBUiKS 



A .cxjmpTiter programme is ■vjritten The assrimption' thqt the 

including all input and output teams temp, of gas at any leve; 
and arSsults presented as graph is higher than thau of 

■between heat and^ temperature, corresponding solid is n 

. ' ■ 3n agreement with, actual 

observation. 




Control of the 
thermal state of 
the blast furnace 


Ore-to-coke ratio rr. - 

< 

• 



Predicting blast 
furnace preformance 
under various 
changes in operating 
va.ria.bles . 

Two: Preparation zone 
in ih ich burden is 
preheated and reduced 
to constite. Elabora- 
tion zone in ■vhich 
direct reduction, 
soln' loss and combu- 
stion, of coke take 
place. 

. r 

J. \ l ' 

111 • i 2 

\ lXi 


Stabilizatiqn of 
Silicon consent 
of pig iron. 


Preparation zone and 
elaboration a)nea 
The dividing horindary 
is an isotherm at 
lOOOoC. 


Blast moist Tire, 

blast temp., msi 

and fuel oil p- rn - 

infection rate ’• 

SolVTi \ iSS 



sini: 

comp 

fern 


Stabilization of 
s ilicon conten t ’ 
of pig iron. 


Preparation '20 ne 
and 

elaboration rone. 


Blast httmidity in i 

■ mat? 
the : 
f-uni 


LongitTidinal dis-. Five: Hearth, . , The 

tribution of temp* Tuyere zone, ■ ' of 

of gas and .solid liquid aane, los 

particles and those fus ion iSeno, stack. sto 

of molar fractions ; 

of CO and C02« 


^5'Co2. -i" 

^^-fO 5*(^C02,'t Hi 0 92 m^.j . j _ /3 A'* £_ . 

Cout.COi.c,i,~¥ d^^)-{7/7(C^^-i-C^^ 
sJSlilirf'' 

M ZXI 2 H H '%S 'Z. ~t 3 " B 'S' 7 

-f 1'-1SA91^) +^eD'3 
=-'7dW4-?^/0 

Z'l z - g^Q 


r- Indices M§ and ^149"^ representing amonr 
of _ neat input in blast furnace per ' 
e^cygen gasified calcu* 
dated from top-gas composition and 
,_from the composition of charge, bias* 
.Irrespectively are closely linked with 
/silicon output of pig iron. 


The blances are made over elaboratic 
33 ne. The third pig. is derived to 
account for limitation arising from 
hea.t .requirement of the preparation 
zone, ' 

The model cannot be used for 
process control. 




ntaneous production; index d^represents the differs 

^ 96 (^A H B) —^2.0 between heat input and the .beat rec 

. -i uirements of solution loss. It is : 

C;i,..ho-n strongly correlated . with the silic( 

~ — .y content of pig iron. The objective 

_ f^y-^ , / ^'7 > to keep Ife as close as possible to. 

P; L’^ ^ ^ h - i'i t (3<-z-tc-o4^\a set value corresponding to the 

©’equations ojre for op^^tion mth _ ^desired ^imlity of pig iron.^This 
er burden. For the operation mth most widely us ed automa.tic con 

lex material similar eqns can be scheme for blast furnaces, 

ulated. 


jidex IS calculated by . The term Ec represents the excess; 

srial and ^ heat ba.lonces worked out in j^eat available in elaboration 20 n* 
elaboration zone. The tryyfer reduce .silica and to superheat 

tion between Ec & Si is j metal and slag to the temp. 

■^corresponding to the silicon cont 
of t^t metal. It is also correlat 
vrith silicon content of pig iron. 


) reacaion rates of indirect reduction 
iron ore by CO and one by Hgj ' 
;s reaction and jiecomposition of lime- 
me were taken into account. 


In this model the temp, in the : 
individual zones is assumed to b 
invariant. 



furnace is higher than that of the corresponding solids at 
that level. The furnace was divided into several temperature 
levels differing by 50°G, and the reactions taking place in 
these zones were taken into account. The heat required by 
the burden passing down the furnace was calculated as a function 
of temperature. The amount of heat available from the gaseous 
phase in cooling from flame temperature to a given temperature 
was taken as proportional to drop in the temperature. The 
resulting curve of heat available versus temperature was nearly 
a straight line. Hodge et al [16~18] presented a scheme to 
predict the changes in the furnace performa.nce provided the 
operating data for the 'standard case' were available. The 
furnace was divided into three regions, the boundaries of which 
were defined by fixing the temperatures of solid phase at the 
boundaries. Three equations were developed, viz., one for 
heat transfer in the shaft, the second for the amount of carbon 
consximed in direct reduction and the third for the amount of 
hydrogen in bosh gas taking part in reduction process. However, 
if the furnace operating conditions were different from those 
of the reference case then the model does not hold good. 
Pokhvisnev et al. [19] divided the furnace into two zoness 
direct reduction and indirect reduction zones. The heat input 

3 

into each zone based on one m at NTP of oxygen removed from 
charge was calculated from top gas composition. The thermal 
state of the furnace was characterized by the difference between 
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the current value of heat input and its reference value. 
Dovgalyuh et al. [20] defined an index 

change in the heat input to indirect reduction zone in the 
time of descent of charge from indirect reduction zone into 
the hearth and was calculated as the difference between the 
heat input into the furnace as a whole and the heat input 
into the lower part of the furnace. Wartmann [21] divided 
the furnace into four zones such as preheating, indirect 
reduction, direct reduction and hearth. The heat and material 
balances in each zone, were combined with the overall material 
balance and the quantity of heat transferred from one zone to 
another, the discharge rate of' slag and the rate of consumption 
of coke were determined. Staib et al. [22-24] pointed out 
the existence of chemical reserve zone and thermal reserve zone. 
The chemica.1 reserve zone is the chemically inert zone in which 
iron exists entirely in the form of wustite. The gaseous phase 
is such that the ratio per cent C02/(per cent CO + per cent CO 2 ) 
is 0.238, a. value set by the equilibriums iron-wust it e-gas. 

The thermal reserve zone is one which extends more or less in 
the shaft where the temperature is periodically constant. The 
chemical reserve zone lies in the lower part of the thermal 
reserve zone. Staib divided the furnace into two parts, viz., 
preparation zone and elaboration zone. The dividing boundary 
vras taken through the chemical reserve zone. The composition 
and temperature of the gas and the solid phases at the dividing 
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boundary were fixed by assuming the temperature level of the 
thermal reserve zone. The thermal state of the "blast furnace 
was characterized by a parameter representing the difference 
between the heat transferred to solids above 1000 and heat 
required for direct reduction per ton of pig iron produced 
and it was strongly correlated with the silicon content of pig 
iron. The objective was to keep close to a set value 
corresponding to the desired quality of pig iron. The model 
is not applicable in cases where the concepts of 'chemical 
reserve zone' and 'thermal reserve zone' are no longer valid. 

The model developed by Daconsine [25] consisted of oxygen 
balance and thermal balance over the elaboration zone. To 
account for the limitations arising from the heat requirements 
of the preparation zone, another equation was also included in 
the model. The model proposed by Vanlangen et al. [26-50] 
was also based on Staib's hypothesis of existence of chemical 
and thermal equilibrium. The thermal state of the blast furnace 
was characterized by a factor B , representing the difference 
between the actual heat supplied and the heat required imder 
standard conditions. It was strongly correlated with the 
silicon content of pig iron. It was calculated by taking heat 
and material balances over the elaboration zone and was 
expressed as the excess of heat available in this region to 
reduce silica and to superheat hot metal and slag above the 
standard conditions, A control scheme was developed using E. 
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as a set point and the blast humidity and blaBttemperature 
as manipulating variables. Flierman et al. [31] divided the 
furnace into five regions such as hearth, tuyere zone, liquid 
zone, fusion zone and the stack. Each zone was treated as a 
stirred tank reactor and the reaction i*ates of indirect 
reduction of iron ore by CO and E^, solution loss reaction and 
decomposition of lime stone were taken into account. Using 
this model, the longitudinal distributions of temperatures of 
solids and gases, and mole fractions of CO and CO 2 were establi- 
shed. The residence time of burden, coke consumption, top gas 
temperature and its composition vrere also calculated using this 
model. Just as Reichardt [14] proposed a graphical model for 
representing heat balances, Rist et al. [32] presented a method 
for calculating both mass and heat balances, by means of a 
diagram known after his name as Rist operating diagram. The 
diagram illustrates most of the chemical characteristics of the 
operation; coke rate, reducing gas composition, gas and charge 
composition at various stages and approach to chemical equili- 
brium. It is also capable of incorporating heat balances as 
constraints on the operating line. The effect of variations in 
various operating parameters like hot blast temperature, oxygen 
enrichment, natural gas injection, prereduction of the burden 
were also studied with the help of Rist diagram. 



Models Based on Differential Heat and Material Balances? The 


mathematical models belonging to this group were presented by 
Koump et al. [33], lahiri et al* [34] and Muchi et al. [35-36] 
and have been sinmnarized in Table 2.4. 

The model developed by Koump et al [33] divides the 
blast furnace into three steady-state chemical reactors and 
one accumulator. The model was developed for Reactor 1 which 
comprised the region of blast furnace from stack line to an 
isothermal surface within the stack where liquid phase begins 
to appear in appreciable quantity. The reactor 1 was treated as 
a steady-state, one dimensional, countercurrent, heterogeneous, 
adiabatic reactor in which components in solid phase are react- 
ing with the gaseous phase. The reaction rates of indirect 
reduction of iron ore by 00 and H 2 , the solution loss reaction, 
the rate of heat and mass transfer between the fluid and solid 
particles , within the solid particles and by bulk flow of gas 
and solid were taken into account. The model consisted of six 
ordinary differential equations pavdlving partial pressure of 
GO and 00 ^ in the gaseous phase, mass flow rates of ore and 
carbon and the temperature of gaseous and solid phases. On the 
basis of the model proposed by Komp et al, , Lahiri et al.[34] 
presented a model for representing both axial and radial distri- 
bution of temperatures of gas and solid particles, mole fractions 
of 00 and 00^ and the fractional reduction of iron ore in the 
stack region of the blast furnace, and considered two reactions - 
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gaseous reduction of iron oxide and gasification of coke. 

Muchi et al. [35-36] developed a model for blast furnace to 
calculate production rate, carbon ratio and the longitudinal 
distribution of process variables in the reactor 1 of Koump, 
Overall reaction rates of indirect reduction Sf iron ore by 
00 and H 2 , solution loss reaction, direct reduction of molten 
wustite by solid coke, decomposition of limestone, reaction 
between coke and steam, water gas shift reaction and reduction 
of silica by solid coke were taken into account. Effects of 
top pressure, diameter of iron ore, volumetric flow rate of 
blast, ratio of steam in]ection, blast temperature, ratio of 
oxygen enrichment and prereduction of iron ore on the producti- 
vity and distribution of process variables in the blast furnace 
under various operating conditions were stadied on the basis of 
this model. 

2.1.2 Dynamic Models ; 

A dynamic model gives the time dependence relationship, 
that IS, the transient behaviour of the process. The develop- 
ment of a control scheme for blast furnace requires a thorough 
study of the dynamic properties of the furnace. The dynamic 
characteristics of a physical system introduces one order of 
complication in a mathematical sense, viz,, if a steady state 
model IS described by a set of algebraic equations, the dynamic 
model will consist of a set of ordinary differential equations 
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and if steady state model is represented by a set of ordinary 
differential equations then the dynamic model will consist of 
a set of partial differential equations. For investigating 
the dynamic characteristics of a blast furnace the method of 
approximation to lumped parameter system is often adopted. A 
large volume element which can be assumed to be completely 
mixed IS selected as a cell instead of infinitesimal element 
to decrease the mathematical complexity. Based on this approach, 
dynamic models were proposed by Fielden et a.1. [37-40], Beer et 
al.[4l], Horio et al. [42] and Tsuchiya et si, [42]. 

Fielden et al. [37-40] divided the blast furnace into 
five regions, namely, stack, upper bosh, lower bosh, tuyere 
and hearth . The furnace was further subdivided into zones of 
one meter height, each acting as a batch reactor. The state of 
each zone was specified in terms of composition and temperature 
of the burden material and gas. A state matrix was defined, 
each row of which described the contents of each zone. It was 
assumed that the gaseous phase and burden material stay in each 
zone for a definite time interval, during which they react and 
move instantaneously to the next zone. An estimate of the top 
gas analysis and temperature was used by the authors to update 
the coefficients of the model. The model was used for automatic 
regulation of blast conditions and charge in order to eliminate 
cast-to-cast fluctuations in the silicon content of pig iron. 
Beer et al. [41] followed a procedure similar to that of 
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Pielden et al. , "by dividing the furnace into n volumes of the 
same size in which physical and chemical changes were deter- 
mined by a system of differential equations obtained from 
heat and mass balances. Horio et al, (see [42]) represented 
the furnace region below melting zone by ' tank-in-series’ model 
and investigated response of lower region of blast furnece to a 
step change of blast conditions. The fluctuations of 'heat 
level' were connected with the magnitude of the coke reserve , 

In the model developed by Tsuchiya et al. (see [42]), the 
furnace was divided into five regions, such as, region 1 in 
which only heat exchange between ascending gas and descending 
burden material takes place, region 2 in vrhich reduction of 
Fe 20 ^ to takes place, region 3 in which reduction of 

Pe^O^ to PeO takes place, region 4 in which gaseous reduction 
of PeO as well as solution loss reaction take place and region 
5 in which combustion of coke takes place. All these regions 
were assumed to act as stirred tank reactors, A steady state 
was assumed with respect to mass balance and with respect to 
heat balance in the gaseous phase. Considering the temperature 
to depend only on the solution loss reaction, the fluctuations 
of solid temperature in region 5 was regarded as a criterion 
of 'heat level and was calculated from observed top conditions. 
This was compared with ihe fluctuations of silicon content of 
tapped metal. Based on a dynamic heat balance, Hodge et al,[43] 
defined three additional variables to supplement the theoretical 
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flame temperature for controlling hearth heat, namely, total 
heat generated , heat available above 2800°F and the ratio of 
heat available above 2800°F to the total heat input. The 
effect of variables like blast temperature, blast humidity, 
oxygen enrichment, fuel oil injection and natural gas injection 
on these three quantities were determined. The authors have 
suggested two more thermal ratios as possible criteria, namely, 
the total heat input divided by the burden charged and the 
heat available above 2800°F divided by the burden charged. It 
was pointed out by Bouman [44] that these thermal ratios do not 
deal with changing reducing requirements in the hearth and 
stack. It IS well known that even with constant blast conditions 
and without apparent changes in burdening, thermal conditions 
are subject to variation, which is indicated by iron analysis 
and its temperature. The change in burdening cannot be detected 
prior to its causing a thermal upset in the hearth. 

Aerodynamic Models ; Bates [45] developed an aerodynamic model 
of the blast furnace relating raw material properties to 
furnace productivity. First, by utilizing classical packed bed 
theory certain relationships like pressure drop, surface area 
etc. were developed, A set of equations relating material 
geometry of the burden materials to the permeability were then 
developed. The material geometry changes inside the furnace 
were related to position inside the furnace, operating conditions 
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and to the different types of material. The furnace was 

divided into five zones, such as, charging zone, indirect 

molten zone 

reduction zone, melting zone,^and turbulent zone and in each 
zone the average values f)r material parameters were used to 
describe that zone, Klempert et al. [46] developed an aero- 
dynamic model for the blast furnace which can calculate 
following process indices; coke retention time in the furnace, 
height of zones of complete and partial melting of materials, 
average lump size for each zone, average size of iron ore 
burden and fine formation in the charge as it descends in the 
region of solid materials and the gas stream temperature up 
the furnace. The furnace was divided into four zones. The 
movable boundaries between zones were calculated by successive 
apnroxima.tion on the basis of material and heat balances for 
one hour of furnace operation. The stack column voidage for 
different zones was calculated with allowance for size degra- 
dation of sinter lumps and their gradual melting as they descend. 

2.2 PROBABILISTIC MODELS: 

In probabilistic models, the technique of statistical 
analysis is used for describing a physical system. The 
probabilistic models are broadly classified as regression models, 
stochastic linear process models and time series models. A 
summary of probabilistic models has been given in Table 2.5. 
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2.2.1 Regression Models ; 

Two approaches have been used in the past to develop 
regression model for the blast furnace, V/hile the first 
approach gives the production rate and coke rate through 
regression models on conourraht values of other variables, the 
second approach gives output variables taking time lag into 
account. 

Por consti-ucting a regression model for operating indices 
(productivity and coke rate), only the characteristics of the 
materials charged and the quality of end products are to be 
considered. With this approach the blast furnace represents 
a 'black box' in which only inputs and outputs are considered 
without knowing the basic mechanism of the process. Such 
models are helpful in production planning and control. However, 
these models do hot give insight to the process and are unsui- 
table for extrapolation and for applying to other similar 
systems. Therefore, these models are, at best, suitable to a 
particular set of operating conditions of a given system and 
outside this range new models may have to be developed. The 
technological parameters of the process that are usually consi- 
dered for regression model are - blast volume, blast temperature 
blast himnidity, consumption of natural gas, oxygen enrichment, 
fuel oil injection, operating intensity, composition of pig 
iron, slag volume, iron content of burden, basicity ratio, 
proportion of sinter in the burden, composition of lime stone, 



coke rate, sinter-to-coke ratio and top pressure etc. Models 
based on multiple regression analysis were developed by Plint 
[47-48], Eudoyarov [49], Starshinov [50] and Kobylyakov etal. [51] 
The models developed by Flint, Kudoyarov and Starshinov were 
linear in parameters whereas Kobylyakov et al. suggested a 
non-linear model. Flint [47-48] expressed coke rate as a 
linear function of 21 variables. Kudoyarov [49] suggested that 
the coke rate and productivity were most affected by the 
composition of burden, and silicon content of pig iron. He 
found that an increase in the blast humidity lowered the 
productivity and increased the coke rate. Starshinov [50] found 
that for an increase in the silicon content coke rate increases 
and productivity decreases. Increase in the operating intensity 
increases both coke rate as well as productivity. Kobylyakov 
et al. [51] developed non-linear relationships between operating 
indices and technological parameters through dispersion and 
correlation analysis and obtained optimum operating parameters. 

The models based on regression analysis with lagged 
variables were developed by Maeda et al. [52], Robbins [53], 
Katsura et al. [54], Spallanzani et al. [55], and Horton [56]. 

Maeda et al, [52] developed an equation relating blast 
humidity with silicon content of pig iron, which was used for 
controlling thermal state of the blast furnace. The blast 
humidity was expressed as a function of amount of direct 
reduction which was calculated from top gas analysis. Robbins[5 



expressed set point for blast temperature as a function of 
blast moisture, wind rate, top gas temperature, carbon monoxide 
content of top gas and heat losses to cooling water. The 
model developed by Es,tsura et al. [54] was based on the fact 
that the difference between the silicon content of pig iron 
and that of preceding tapping is directly proportional to the 
heat storage in the lower part of the furnace. This difference 
of the silicon content was related to the variables like blast 
temperature, humidity, flow rate; amount of direct and indirect 
reduction; amount of heterogeneous reaction and permeability 
in the furnace. Spallanzani et al. [55] assumed a hypothesis 
similar to that of Katsura et al. and predicted the silicon 
content and hot metal temperature of next cast based on variables 
measured at 9 and 5 hours before the expected casting time. 

While Robbins had completely ignored the time lags (between the 
charge and the cast and between the blast and the cast), Katsura 
et al. have considered the time lag between charge and the 
cast and Spallanzani et al. have considered both the time lags 
and constructed the model by delaying: certain measurements. 
Norton [56] has developed a model in which a.n output variable 
at any time is expressed as a linear function of the input and 
output variables at previous sampling times and a noise. The 
equation was linear but the variables which had non-linear 
relationship with the measured variables were separately 
estimated and used in the linear regression. The linear models 
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were developed for silica reduction rate and production rate 
by linearization about the normal operating conditions. The 
production rate model was based on material balance and the 
silicon reduction rate model was based on the heat ba,lance 
above 1000°G. 

2.2.2 Linear Stochastic Process Models ;.- 

Another group of the probabilistic models is the linear 
stochastic process model. The input and output variables of 
the blast furnace are subject to random perturbations that can 
be modelled by probabilistic laws and hence a linear stochastic 
process model can be used to relate input and output variables. 

A steady-state is first carefully established, A 
specified perturbation (impulse, step or sinusoid) can be given 
to an input variable and the response of the system on one or 
more output variables may be measured. They are used in esti- 
mating the system pa,rameters. This is an ’ instrimnentation 
problem’, luckers et al. [57-58] used binary perturbations, 
ma.inly square waves to estimate coefficients of various 
tra,nsfer functions, whose forms were chosen from previous 
experience. Models relating blast temperature, blast humidity, 
oil injection rate and coke-to-sinter ratio to the silicon 

content of the hot metal and the narameter B were fitted. 

" c 

After logging the furnace response at a sampling interval of 
20 minutes the cross-correlation function between input and 



output was computed. The plant model was then fed with a 
waveform representing the autocorrelation function of the 
perturbation input and the coefficients of the transfer function 
were adjusted until a satisfactory fit was obtained between 
the model output and the recorded cro ss~correlation function. 

The ¥iener-Hopf equation was thereby used in which knowing the 
autocorrelation function of the input , the impulse response 
was adjusted until best input-output cross-correlation function 
was obtained. It was found that silicon response to blast 
temperature variations and changes in the flow of fuel oil is 
very slow. The transfer function between blast humidity and 
silicon content had shorter time constant but its form was’invers< 
response type* [58], The silicon response to the coke-to-sinter 
ratio included a dead time in addition to the time constants. 
Staib et al. [59] have also studied the effect of step response 
of fuel oil rate on hot metal silicon content. Rebeko et al.[60] 
also used Wiener-Hopf equation to determine the impulse response 
function knowing the input-output cross-correlation function 
and the input autocorrelation function. The mathematical model 
was based on calculation of available heat in the direct 
reduction zone at a series of successive moments of time after 
a step wise change in ore-to-coke ratio. The system transfer 
function was obtained from the impulse response function. The 
transfer functions between the variables like ore-to-coke ratio, 
iron content in sinter, oxygen content of the blast, proportion 



Sf tn^'fOiral gas in the blast and silicon content of pig iron 
were determined. • 

The classical methods of estimating transfer functions 
based on deterministic perturbations, however, have not always 
been successful because of the following reasons? (i) These 
tests require carefully controlled furnace conditions for 
several days per test and create considerable operating incon- 
veniences. Similarly the sequence of pseudo-random binary 
signals would each require several days of data logging, during 
which a few missing or unreliable points or necessary control 
action would spoil the experiment, (ii) The blast furnace is 
subjected to a large number of disturbances which cannot be 
taken into account in such a model. (iii) It is difficult to 
decide on the amplitude of variations to be imposed so that it 
should take into account, on one hand, the size of the distur- 
bance and on the other hand the upper limit imposed by the need 
to remain close to the normal operating. conditions. Some of 
these limitations are overcome in time series models. 

2.2.3 Time Series Models ; 

Time series models were developed by Rylov [61], 

Rake [56], Rebeko et al. [60], Gastore et al. [62] and Phadke 
et al. [2]. 

Rylov [6l] applied the methods of statistical dynamics 
to study the relationships between rate of descent of charge 
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material with, the input quantities - blast volume, temperature, 
humidity, and natural gas consumption. On the basis of cross- 
correlation coefficient it. was shown that the relationship 
between the rate of descent of charge and blast volume had 
an extremum. A system of - extremal regulation for the rate of 
descent of charge by changing the blast volume was proposed. 

Rake (see [56]) determined the transfer function at a given 
frequency between the heat index analogous to the parameter 
and the ■ silicon content of pig iron from the power spectrum 
of the input and the cross-spectruun between the input and output 
These spectra were calculated by taking Fourier transformation 
of autocorrelation function and input-output cross-correlation 
function. Rebeko et al. [60] used the methods of statistical 
dynamics and the theory of random functions for predicting 
thermal state of the blast furnace, particularly silicon 
content in the next cast, fromthe informe.tion on the current 
operation of the furnace. On the basis of primary information 
for a definite time interval, the autocorrelation and cross- 
correlation functions between input and output variables were 
determined. The position of cross-correlation maxima was 
determined and the coefficients of multiple regression equation, 
between the parameters taken with appropriate time displacements 
were determined. Castore et al. [62] have developed a feed- 
back control scheme for hot metal silicon content using blast 
humidity, blast temperature and oresto-coke ratio as 
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ina.nipulating variables. The input variables were subjected 
to programed stepwise variations and transfer functions 
between these independent variables and silicon content were 
developed. The disturbance at the casting frequency was 
identified by means of statistical analysis of time series 
of finite length collected during the periods of normal, 
conditions. In the feed back loop, the output of the ’plant 
model’ after sampling and delaying, and the disturbance which 
represented the variable component of silicon content in 
open loop, were fed to the controller which determined the 
correction on the control variable. Phadke et al. [2] have 
developed a multivariate time series model for the blast furnace, 
taking into account two input variables (blast flow rate and 
ore-to-coke ratio) and three output variables (sulphur content, 
silicon content and temperature of molten metal). Using this 
multivariate model (described in detail in Chapter 4), they 
have also derived various transfer functions, viz., for the 
plant, for the feedback system, for the blast furnace noise and 
for the feed back noise. 

Because of the advantages of time series modelling over 
other typs of models for the blast furnace process, it is 
proposed to adopt time series modelling in this study. In 
particular, Phadke’ s model is adapted and used. 


+++++++++ 



CHAPTER 3 


UNIVARIATE TIME SERIES ANALYSIS 


3.1 INTRODUCTIONS 

A time series is defined as a set of observations made 
sequentially in time. If the set is continuous, the time 
series is said to be continuous. If the data are given for 
discrete values of time as sampled or quantized data, the time 
series is referred to as a discrete time series. If the future 
values of a time series are exactly determined by a mathematical 
function using the present and the past values, the time series 
is said to be deterministic. If the future values can be 
described only in terms of a probability distribution, the time 
series issaid to be a non-deterministic or a stochastic time 
series. 

A time series is said to be strictly stationary if the 
joint probability density function of the families of the 

variable [z(t^), z(t 2 ), ,z(t^)] and [z(t^_^j^), z(t 2 ^j_), 

z(tj^_^i] is the same for all i for any set of t^, . ,tj^; 

ie. , stationarity is achieved if the distribution is invariant 
with respect to translation in time. This implies stationarity 
of moments of all orders of the time series. Eirst order 
stationarity implies that mean is constant and second tErrder 
stationarity implies that the auto covariance at I9g k, viz. , 



E[z(t) ’zCt+k) ] , where E stands for the expectation, is a 
function of k only. This also implies that the variance of 
z(t) is constant, G-enerally a second order stationarity is 
assimned. 

A time series can generally be considered to consist 
of deterministic components that include trend and cyclicity; 
a non-random component which includes peristence^ and a pure 
random component. Trend refers to long-term behaviour of the 
time series and it may indicate a mono touQU a change in the main 
value - which can be expressed as a linear, polynomial or 
exponential function of time, and/or jumps in the mean value. 

1 plot of data may indicate whether a trend is present. 

Cyclicity refers to the periodic components that repeat them- 
selves at definite intervals. Persistence refers to the linkage 
between the value at 'a given time with earlier values and may 
be due to internal or external dependence, A time series may 
be represented in terms of its several components as follows: 

z(t) = ftrend^"^^ ^cycle^^^ f^j^(z(t-l) ,z(t-2) z(t-p)) 

+ f^_^(e(t) ,e(t-l) e(t-q)) + eit) 

( 3 . 1 ) 

in which e(t) represents the pure random component of the time 
series and f^^^ and f^ represent autoregressive and moving 
average components, representing respectively the internal and 
external dependence in the time series. The random component 
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is the residual series obtained after subtracting trend, 
cycle, and autoregressive and moving average components from 
the original time series. Very often some components may not 
be present and the resulting time series is simpler. 

5.2 TIME SERIES MODELS: 

An important class of stochastic models for describing 
a time series is the so called stationary models. A time series 
is said to be stationary if the generating mechanism of the 
process is invariant with time. This implies that the deter- 
ministic component of the series is independent of time and 
the parameters of the model are also independent of time. 

5.2.1 Sum of the Harmonic Process Model ; 

In order to identify deterministic components of a 
time series, it can be decomposed into a number of sinusoids 
of varying frequencies and amplitudes. When the data are 
finite and discrete with N=T/ ^t samples, where T is the time 
horizon and At is the time interval, z(t) can be approximated 
by a finite Eourier series that passes through all sample 
points, viz., 

M 

z(t) = cj^ sin(w^t + ) + x(t) (3.2) 

i=l 

where c^ represents amplitude, is the phase angle and w^^ 
is the angular frequency given by 


2 71 f^ = 2 IT i/T 


(3.5) 



in which f=l/T is the fundamental frequency and i represents 
the ith harmonic. If the shape of the time series is. that of a 
sine curve then the entire fluctuations of. the .jseriea would 
be contained in one hamonic. When, the shape of- function is 
not sinusoidal, then atleast two and atmost N/2 -harmonics are 
required to approximsite a function through all sample points. 

3 . 2.2 Autoregressive Model ; 

This is a linear model for persistence within the time 
series, in which the current value of the process is expressed 
as a finite, linear aggregate of previous values of the process 
and a random component. It is given by . .. 

x(t) = 9^x(t-l) + 92^(t-2)+. . . .+9pX(t-p)+ e(t) 

■ •(■ 3 . 4 ) 

where x(t) represents the deviation from the mean p of the 
original time series x(t) and e(t) is a random component of 
mean-, zero. This is an autoregressive process of order p or 
AR(p) process. The backward shift operator B is defined by 

Bx(t) = x(t-l) and B^x(t) = x(t~k) ( 3 . 5 ) 

Then Eqn, 3.4 can be written as 

9(B) x(t) = £ (t) ■ ( 3 . 6) 

where 9(B) is the autoregressive operator of order p defined by 

9(B) = 1-9^B - -9pBP ( 3 . 7 ) 

Let represents the variance of the white noise process e(t), 
then the model contains p+2 unknown parameters, 9^,92,....cpp 



and which ha.ve to be estimated from the data, 
e 

3.2.3 Moving Average Model ; 

This is’ ja'‘niQd.e3jdind±cating external correlation between 
the, time series x(t) and a pure random time series e(t). In 
this model, the current value., of the process is expressed as 
a linear function of the present and the previous values of the 
random component e(t), and is given by 


x(t) = e(t) - O^e(t-l) - 02 e(t- 2 )-, . .-&g^e (t-q) (3.8) 


This is a moving average process of order q or MA(q) process. 
Equation 3.8 can also be written in terms of backward shift 
operator, B, as 

x(t) = 9(B) e(t) (3.9) 

( 3 . 10 ) 

is the moving average operator of order q. The" model contains 


where ©(B).^.^. 1 - © B - ©pB^- 


q+2 unknown parameters p; ©., , ©^,...., © and which have to 

X ^ q c 

be estimated from data. 


3.2,4" Autoregressive - Moving Average Model ; 

To achieve greater flexibility in fitting a time series 
model, it is sometimes advantageous to include both the 
autoregressive and moving average components in the model. 

This leads to a mixed or autoregressive-moving average(ilBJlA) 
model of order (p,q); 

x(t) = 9 ^ x(t-l) + (p 2 ^(t- 2 ) + , . .+(ppX(t-p) + £(t)- ©^e(t-l) - 

©qe('t-l) (3.11) 




■which can also he written as 


cp(B) x(t) = 9-(B) e(t) (3.12) 

It contains p+q+2 •unknown parameters that have to he estimated 
from the data. 

3.2.5 Autoregressive-Integrated -Moving Average Model : 

The above models are stationary ■with respect to mean 
and standard deviation and they also assume covariance statiou 
arity, A more general class of models is one which is non- 
stationary and whose stationarity can he achieved hy repeated 
differencing of the original series. Let 

w(t) = x(t) (3.13) 

where is the backward shift operator defined as 

^x(t) = x(t) - x(t-l) (3.14) 

and d is the order of differencing. A stationary AEMA(p,q) 
model is then fitted to the differenced series w(t). The 
original series x(t) is said to he autoregressive-integrated- 
moving average (ARIMA) model of order (pyd,q) where d indicates 
the order of differencing required to render the original 
process stationary. 

In general, the time series, can he expressed hy the 
general ARIMA model in which p and q are not greater than 3 
and d is not greater than 2, 



3.3 IIESAR SYSTEMS; 

A system, whose behaviour can be described by a linear 
differential or difference equation is said to be linear. 

3.3.1 Linear Differential Equation ; 

A continuous linear process with input e(t) and output 
x(t) ca,n be represented by a linear differential equation 
with constant coefficients, as follows; 

[a +a D + a D^+,..,+a D^l x(t) = 

0 1 2 p ' 

[ b^+b^D+b2D2+ .... +b^D'l ] e ( t ) (3.15) 

where a's and b’s are constants and D is the differential 
operator. A system characterized by a differential equation 
is dynamical, because the value of the output x(t) at any 
instant t depends not only on the value of sCt) at that instant 
but also on the values of the derivatives of £(t) and x(t) 
at that instant which, in turn, are dependent on the values of 
x(t) and e(t) at other instants. A system described by Eqn.3.15 
is time invariant if the differential equation relating the 
input and output signals has constant coefficients. 

3.3.2 Linear Difference Equation ; 

Just as in the case of a continuous system, many 
properties of a discrete system are exhibited in the difference 
equation relating the input and output signals of the system. 

A discrete linear system with input e(t) and output x(t) may 
be represented by a linear difference equation with constant 




coefficients, as follows? 


c^x(t) - c^i(t-l) - G2x{t-2) -CpX(t-p) = 

v^e(t) - v^e(t-l) - ^ 26 ( 1 - 2 ) . . . ,v^e(t-q) (3.16) 

where c's and v's are constants. 

A discrete system described by Eqn. 3.16 is dynamical 
because the value of x(t) depends not only on the values of 
£(t) but also on the values of x(t-l), x(t-2) , . . . . ,E(t-l) , 
e (t--2) , . . . . . The system described by Eqn. 3.16 is also time 
invariant if the coefficients of the difference equation are 
all constants. The general AEMA model given by Eqn. 3.11 
is exactly similar to the linear difference equation given by 
Eqn. 3 . 16 . It can be seen that there is close resemblance 
between the linear differential equation with constant coeffi- 
cients and the ABMIi model. 

3.3.3 Stability ; 

A system is said to be stable [ 63]» if the magnitude 
of the output is bounded at all times when the magnitude of 
the input is bounded at all times. A system which generates 
the time integral of its input need not be stable since the 
magnitude of the output can be unbounded even when the magnitude 
of the input is bounded. A discrete system whose output is 
the forward or backward difference of its input is always 
stable. Since a pure moving average process, described by 
Eqn. 3.8 can be considered as an output from a linear system 
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whose input is a white noise, it is always stable. 

Stability of OOphtinuous System ; The characteristic equation 
of the differential equation 3.15 is as follows; 


a s^ + . . . . + a, s + a^ = 0 (3.17) 

p 1 0 

The system is said to be stable if the roots of the character- 
istic equation have negative real parts. 


Stability of Discrete System ; let z be a shift operator with 
the property 

z^^xCt) = x(t-k) (3.18) 

Then the difference equation 3.16 may be written as 


(o^-o^ x(t) 


(v +V, z“^ + V„z“^+. +V^z“^)£(t) 

0 1 2 q 


(3.19) 


thus 


V ' '' T 

Xit) = 


v_+v^ z“^ +V2z"^+. . .+Vg^z"^ 


c -C^z ^-c„z -....-c z ^ 
o 1 2 p 


(t) 


( 3 . 20 ) 


The characteristic equation of a discrete system is obtained 
by equating the denominator of Eqn. 3.20 to zero. Thus 

...-c^ =0 (3.21) 


C Z^ ~ CSt 
0 1 


The system is stable if the roots of the characteristic 
equation lie inside the unit circle in z plane. 
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3.4 ANiiLYTICAI PROCEDURES: 

The periodic and persistence components of the time 
series are identified hydcorrelation analysis and spectral 
analysis. The autocorrelation function, partial autocorrelation 
function and spectral density function are important tools in 
identifying the periodic and persistence components. These are 
helpful not only to identify the form of the model but also to 
obtain approximate estimates of the parameters. 

3.4.1 Autocorrelation Function ; 

The autocorrelation function determines the linear 
dependence among successive values of a series that are a given 
lag apart. The auto covariance function between two values z(t) 
and z(t+k) of a time series that are tk lags apart is given by 


Y(k) = cov[z(t ) .z(t+k) ] =E[ (z(t )-p(t ) ) (z(t+k)-ij.(t+k) ) ] 

( 3 . 22 ) 

The autocorrelation function, at lag k is given by 

£ _ pQ'^r z(t ) «z(t+k) 1 _ Ef (z(t)-u(t) ) (z(t+k)~u(t-fk) ) 1 

[varz(t) •varz(t+k)]^ ^E[ ( z (t ) -p ( t ) ) ^ [ E[ ( z ( t+k ) -p ( t+k) ) ^ ] 


(3.23) 


The plot of the autocorrelation function with lag is known as 
a correlogram. 


Estimation of Autocorrelation Eunction ; In loractice ,from a 
finite sample of the tim.e series one can only obtain the estimate 
r^, of the autocorrelation X-^t where 
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'k 


in which 




R 

0 


(3.24) 


k 




‘ (z(t)-(i ) ( 2 (t+k)- ) 


is the estimate of autocovariance For computational 

purposes, the autocorrelation function may he calculated from 

F-k 


F-k i=l' 


F-k F-k 

>^z(i) z(i+k) n_ [( ^ z(l))-(y’a(l+k))] 


irST 


i=l 


(var z(i)*var z(i+k))' 


F-k 


1= 



(3.25) 


where var z(i) = ^ z(i) - [ / 

i=l 


z(i) ]' 


Standard Error ofBstimated Autocorrelation Function ; The 
estimated value of the autocorrelation function differs from 
the theoretical value because of sampling errors and finite 
sample size. Hence, it is important to have some indication 
of how far.rin estimated value may differ from the theoretical 
va,lue. In order to test that the autocorrelations are not 
significantly different from zero after some log q., Bartlett[64] 
has suggested the following equation for the standard error of 
estimated autocorrelations, 

“V“ [1+2 (r^ + r^+ ^ ^ (3.26) 

^ ]^2 ± ^ 4 . 

so that all autocorrelations beyond + 1.96 times the standard 
error may be considered significant at 95 per cent confidence 


level. 



3.4.2. Partial Autocorrelation Function? 


The partial autocorrelation function is useful in 
identifying particularly the order of an autoregressive 
process. The partial autocorrelation coefficients are related 
to the autocorrelation coefficients by Yule-¥alker equations 
given by [see 65,66] 
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(3.27) 


If the order of the process, p, and the autocorrelation 
functions, known, the above system of p 

equations with p unknowns (p2_> 92 »**'*» 9p solved. In 


practice, the true values of p and are unknoTAra. Let p=l^ 

Using Yule-¥alker equations and using the estimated value r^, 
for one gets r^ = 9-,, , where is the estimate of 9^. If 

9^ is significantly different from zero, it can be concluded that 
the process is atleast of order one. To see whether the process 
is of order two or greater, the Yule-Walker equations are solved 
for p=2, i. e. , 


^1 = 


9i + 92 ^1 


9 i 


9 -: 


( 3 . 28 ) 


in which r_ stands for the estimate of theoretical autocorrelatio 
k 

coefficient, 





If the resulting estimate of cp2 differs significantly from 

zero, it can he concluded that the process is atleast of order 

2. This procedure is repeated successively for larger values 

of p. If the true order of model is fhen, when the 

system is solved for p=p.„„ +1, the value of will not 

® ^true+1 

he significantly different from zero since it is an estimate 

of <p , which is zero. Denoting hy cp. . , the Talue of cp. 

' ^true+l ^ ^ 

obtained hy the solution for p=i, 9^^ are referred to as the 

estimated partial autocorrelation coefficients of the process. 


If the order of autoregression is then 


= 0 for 


i > P 


•‘^true 


(3.29) 


The partial autocorrelation coefficients are calculated using 
Ec[n, 3.30 


'Pii= i 


where 9- • 

-L ^ 




i-1 

, . r. 


i= 2 , 3 ,. . .. 


( 3 . 30 ) 


‘P\-1,0 " ^ii ?i-l,j-l 3-1,2,.. .1-1 


Standard Error of Estimated Partial Autocorrelation Goefficients s 
Since the partial autocorrela^tion coefficients are sample 
statistics and therefore subject to sampling error, a test is 
needed to decide when “9'^^ is indistinguishable from zero in a 
statistical sense. Quenouille [6 7 ] has shown that on the 
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hypothesis that the process is autoregressive of order p, 
the estimated values of order p and higher are approximately 
ncrmally distributed with standard deviation 

- 1/VN k p+1 (3.51) 

Thus it is inferred, that p = Pj^. at 15 per cent confidence 
level if 9 is small compared to + 1.96/V^N. 

Pk 

3.4.3 Snectral. Density Function ; 

Spectral analysis is useful in identifying periodic 
and persistence components of the time series. For discrete 
data the frequency domain representation of z(t) is a discrete 
line spectrum,. A complex periodic data consist of a number of 
sinusoids with amplitudes c^ , and phases [see Eqn,. 3.2], 

The plot of amplitude versus frequency is called a periodogram. 

The spectral density function is a normalized periodogram and 
it shows how variance of z(t) is distributed over the frequencies. 
It is obtained by taking Fourier transformation of. the auto- 
covarance function. If M represents the maximum lag (usually 
10 per cent, of the number of data points) upto which auto- 
covarance, is computed, then the sample spectral densitj^ 

function is estimated by 

M-1 

G^ = G(f ) = I [R^+2 ^ R^ cos + (-1)^' Rj^] (3.52) 

where f is the frequency given by 
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f — 2M zi t * i — 0, 1, 2j . . . ,M (3.33) 

The cut-off frequency, f = l/(2*At) is the hipest frequency 
cycle that can be identified from a given discrete series 
with time interval ZAt. The spectral density function G-Cf) 
has a negative exponential distribution and so the sample value 
G(f) may be widely different from the true population value. 

This is because , while the sample autocovariance function is a 
good estimate of the point value, the simultaneous estimates 
for all points are not good estimates. It is difficult to 
identify the components of the time series from the raw spectrum. 
Hence, it is necessary to smoother the raw spectrum. There 
are several procedures available for this purpose. One way 
is to compute the spectrum and then perform the smoothing 
operation. Another way is to multiply the autocovarance 
function by a weighting function, a.nd then compute the spectrum. 
In this study, a method suggested by Blaclonan and Tukey[see 68] 
is used, viz,, smooth spectra are given by 

G = 0.54 G^ + 0.46 G-, 
o o 1 

Gj^ = 0.54 Gj^ + 0.46 (3.34) 

G, = 0.23 G, , + 0.54 £ + 0 . 23 % 1< k <M 

k k-1 k k+1 

3.5 STEPS IN PITTING A UNIVARIATE TIME SERIES MOEEL; 

There are three steps in fitting a time series model, 

viz. , 
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(i) Identifications Using the data the order of 
differencing and the form of the model are determined. 

(ii) Estimation: Initial (approximate) parameter 
estimates are obtained from autocorrelation function and are 
used as starting values. The parameters of the identified 
model are then estimated using more refined procedures. 

(iii) Validation of the Model: The residuals from the 
fitted model are subjected to diagnostic checking to test the 
adequacy of the model. If the residuals are purely random, 
no lack of fit is indicated, and the model is considered to 
represent the physical system. If any inadequacy is found, 
the iterative cycle of identification, estimation and diagnostic 
checking is repeated until a suitable representation is found. 


5.5.1 Data Used for Study ; 

Pig iron is produced in blast furnace by reducing iron 
ore. The reducing agents are hydrogen, carbon monoxide and 
carbon. The reducing gas is produced by combustion of coke by 
blast air at about lOOO^C and decomposition of steam coming 
with blast air. The reactions are exothermic and produce large 
amount of heat in the lower part of the furnace. The reduction 
of iron ore by CO and is called indirect reduction and that 
by solid carbon is called direct reduction. If the reduction 
is more than noimal by either of the process, it causes 


imbalance in the operation, as a result of which the temperature: 


of the hearth increases. Since the solubility of. silicdhltin the 

CENh'.AL liBRARY 

T’^nfST** A 


Mmm Ho* 
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molten metal increases with the temperature silicon composition 
of the hot metal would increase. If the silica in the burden 
IS increased, other conditions remaining uniform, silica content 
of slag will increase whereby corresponding increase inthe 
silicon content of molten metal will be noticed. 

The two Tariables used for controlling hearth heat and 
hence the silicon content of hot metal are blast humidity and 
sinter-to-coke ratio. The effect of these two variables on 
silicon content are as follows: 

Blast Humidity: It immediately reduces the temperature of the 
hearth due to decomposition of water vapour and correspondingly 
there will be a reduction in silicon content of the hot metal. 
However, increase in humidity increases the reducing power of 
the gases due to hydrogen andhence indirect .reduction is favoured. 
Thus more carbon is available at the hearth for direct reduction 
which increases the temperature of the molten metal and hence 
its silicon content. 

Sint er-to -Coke Ratio; When sinter-to-coke ratio is increased, 
the supply of coke per unit of sinter is decreased. The availa- 
bility of carbon for direct reduction, in the lower region of 
blast furnace is reduced and hence hearth temperature is lowered. 
This will reduce the silicon content of hot metal. The physical 
properties of slag will also be affected. The long term supply 
of the reduced quantity of coke will depend upon the slag 
properties. The decrease in the coke may partly be compensated 
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by increasing the blast flow rate. 

The input variables have long term and short term 
effects on the operation of a blast furnace. However, their 
effects are insignificant if the magnitude of variation is 
limited because of the high capacity of the system. 

The data on the following input and output variables 
have been collected on Blast Furnace Ho.l of capacity 2000 tons 
per day of Bokaro Steel Limited, Bokaro Steel City, 

Input Variables; 

(i) Sint er-to -coke ratio ■' 

(ii) Blast flow rate 

(iii) Blast humidity 
Output Variables; 

(i)'Hot metal temperature 

(ii) Silicon content of hot metal 

(iii) Sulphur content of hot metal 
The metal was cast approximately 10 to 11 times per day with 
the cast period varying from two hours to three hours. Data 
have, been collected for 350 consecutive casts comprising the 
operation of 36 days. The input varie.bles were recorded 
continuously and sampled at one hour intervals. The time series 
■ for input variables were ' generated by averaging the values 
between successive casts. The observations on output variables 
were available only at the end of each cast period. Though 
the cast periods were not equal, for simplification in the 
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analysis, it has been assmaed that all periods are equal and 
the cast period is considered as a 'time step*. The data have 
been reported in Appendix G and are shown in Figures 5.1 to 
5.6. It should be noted that values of hot metal temperature 
were not available at all the casts. These points have not 
been included in Figure 5.4. 

5.5.2 Preliminary Analysis of the Data : 

All six variables, viz., sinter-to-coke ratio; blast 
humidity; blast flow rate; temperature, silicon content and 
sulphur content of hot metal have been named A,B,C,D,E and F 
respectively. It is seen from Figures 5.1 to 5.6 that there 
are two jumps in the series C and so it is nonstationary. 

Hence, it was divided into three parts, namely, and 

C^, each of which was treated as stationary. All three parts 
ha,ve different means and standard deviations. All other series 
are apparently stationary. The means and standard deviations 
of the variables are given in Table 5,1. 

Out of 550 casts, the data on hot metal temperature 
were available only for 270 casts. Hence, the model was 
fitted to hot metal temperature with a. sample size of 270 and 
then using the fitted model, the missing values were estimated 
and used in parameter estimation. This procedure was repeated 
until there was no change in parameter estimates in two success- 
ive iterations. Since the order of magnitude of mean and 
standard deviation is different for all series, the series were 
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Table 3.1 

MEAN AND STANDARD DEVIATION OE DATA SERIES 


Mean 


Standard 

deviation 


2.458 


7.161 X 10 


-2 


iSeries 


No . of 
observations 


Sinter-to-coke 
ratio, A 

350 

Blast humidity, B 

350 

Blast flow 
rate, 0^ 

43 

Blast flow 
rate, C 2 

180 

Blast flow 
rate, 

127 

Hot metal 
t emperature,D 

270 

Silicon content 
of hot metal, E 

350 

Sulphur consent 
of hot metal, E 

350 


45.54, gm/m^ 

3.841 

3111.0, m^/min 

18.38 

3411.0, m^/min 

27.53 

3237. C, mVmin 

18.57 

1455. 5, °C 

28.05 

1.404, per cent 

0.0235 

3. 9xl0~^,per 

8.29x10 


cent 


standardized using '\x the mean and "a the standard deviation 
of the raw series, viz., 

y(t) = [z(t ) - VI ]/o 

so that all transformed series have zero mean and unit standard 
deviation. The standardized series are then used for further 
analysis. 

3.5.3 Identification ; 

The aim of identification is to identify the presence 
of trend, persistence and periodic components in the time 
series. 

Trend ; The trend is identified by visual observation on plott- 
ing the raw data. If a trend is present, it can be estimated 
by regression analysis. A linear trend can be removed by 
taking the first difference of the series. The estimated trend 
component is subtracted from the series to give a trend free 
series which can be used for further processing. 

Periodicity ; The presence of cycles can be identified from the 
plot of the series and also from the autocorrelation function 
and spectral density function of tihe series. The presence 
of cycles is indicated by the occurrence of periodicity in the 
autocorrelation function and spikes at appropriate frequencies 
in the spectral density function. Once the frequencies are 
identified, the amplitudes and phase angles are estimated by 
harmonic coefficients or by regression analysis [69 ]. The 
cycles, after being identified, are subtracted from the trend 



free series tc give a series free from cyclicity. 

Persistence ; It refers to linkage ‘between the values at a 
given time with earlier values and may "be due to internal and/ 
or external dependence. The trend free and cycle free series 
IS analysed for the presence of persistence using autocorrelation 
function, partial autocorrelation function and spectral density 
function. The general behaviour of these will reveal whether 
the process is autoregressive (internal dependence), moving 
average (external dependence) or mixed. The autocorrelation 
function of an autoregressive process of order p tails off and 
its partial autocorrelation function has a cut-off after lag p. 
Conversely, the autocorrelation function of a moving average 
process of order qfhas a cut -off after lag q, while its auto- 
correlation function decays gradually. If both autocorrelation 
function and partial autocorrelation function tail off, a mixed 
process is indicated. 

Identification of the Plant Data ; The autocorrelation functions 
for all series vrere calculated usingEqn. 3.25. The partial 
autocorrelation functions and the standard errors were calculated 
usingEqn. 3.30 and Eqn. 3.31 respectively. The autocorrelation 
and partia,! autocorrelation functions, along mth the 95 per 
cent confidence level for the latter have been plotted in 
Figures 3.7 to 3.12. 

The raw spectral density function and smooth spectral dens 
function were calculated using Eqn. 3.32 and Eqn. 3.34. The 
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snooth. spectra, for all series, have been plotted in Figures 
3.13 to 3.18. They indicate no periodic components in any 
series. Based on the characteristics of autocorrelation 
function and partial autocorrelation function, a tentative 
identification of the models is given in Table 3.2. For most 
of the series more than one model was identified and the 
model with best f'lt was selected later at the stage of diagnostic 
checking. 


3.5.4 Estimation of Parameters: 

Having tentatively identified one or more models for 
a time series the next step is to obtain on the basis of some 
criterion, the best estimates of the parameters of the models. 
To start with, approximate estimates of the parameters are 
required. These are obtained from the autocorrelation coeffi- 
cients of the process as followss 


Autoregressive Process. AR(g) i For an ARijp) process, the initial 
estimates are obtained by solving the Yule-Walker equations 
(see [ 65 ] and [66])- 




•p-1 

f> 

■p-2 


r ^ r . 

p-1 p-2 




^1 

A 

92 


% 


(3.35) 


The initial estimate of the residual variance is obtained from 


- 


(3.36) 



Spectra G (f) Spectra G(fJ 



-Power spectra of sinter- to-ccke ratio senes A 
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where, is the variance of the original series. 


Moving .tiverage, MA(q) Process ; For a MA(q) process the initia 


estimates for q unknowns 
the following equations j 




, Q are obtained from 
1 


A /s ^ y\ 

— , n + . • . 0 +© T © 

r = — ^ OrL-i , k=l,2, q (3.37) 


(1+ 0'^ + ©; + + 

X ^ q 


The initial estimate of the residual variance is obtained from 

R 


£ 


0. 


(l + + ©2+ • • • • 


( 3 . 38 ) 


Mixed. ARMA(o.q) Process ; In the general case, the calculation 
of initial estimates of an ARI'([n(p,q) process is based on the 


first Cp+q+l) autocovariances R [3=0,1, .... (p+q) ] of stationar 

J 


series z(t). The autoregressive parameters 9^, cp2s • • • • » 9 p 


estimated from the autocovariances R_ „ , ,R , , R 


q-p+l» 


q+l» 


by solving the following set of equations: 


q+p 


R 


/"V 


R 




q+1 'f’l q "^2 q 
'q+2 ^ ‘Pl\+1 ^ ^2^ 


. . , + cp R^ , T 
q-p+i 

P \-p+2 


+ 9, R 


( 3 . 39 ) 


q+p ^1 q+p-1 ^2 q+p -2 ^p q 

Using the estimates obtained from Bqns !|. 39 » the first (q+l) 

autocovariances R' , (3=0,1, ... .q) of the derived series 

3 

Z*(t) = z(t) -9 z(t-l)y...- 9 z(t-p) 

ir 

are calculated from 
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) iL feo ®fO+i-k| 


R’ 

3 


_ \ 


P > ° (9^0=-^) 


00 


) 


(3.40) 


R 


p = 0 


XV 


where ;) = 0,l,,o.,q a,nd sta.nds for initial estimate of 

ith AR parameter. The initial estimates of the moving average 
parameters are obtained by a method suggested by Wilson [70] 
[65]' using Newton-Raphson algorithm. Let 

T 

r ” ( t ) P-j y « * . . ) 

_ 0 1 q 


where ; and 9 = -t /r 

o e’ 3 30 


3 = 1,2,... ,q (3.41) 

Then, if is the estimate of t obtained at ith itero.tion, 

the new values at rhe (i+l) iteration are obtained from 


X =1 


j 1 ^ ^,1 


(3.42) 


q-3 


and 


T = 


f* = 


■p 

9 -A T J • 

f ' )'^ 

...iq; , 

f' : 
3 


\ 

i=:0 

T 

X 

T ^ -E' 
1+3 3 



T 

X 0*9 

'T 

T 


X 

T- 


9 9 X 


0 


•** q-1 

q 


0 

1 

2* ’ ’ 
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% 



0 


0 

z 

T 9 9 

T 


1 

T 2 9 9 0 

. . . . tq 


0 

f. . ' . 

•• q-1 

== 

x^ 

T 

. . . ,0 

0 


0 

0 

X 9 9 9 9 

9 9 1^ o 


2 





0 

q-2 


'^k 

0 

0 

0 


0 

0 

0 

0 


(3.43) 
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The starting values are ^ , '^1 ~ '^2 ~ ..-.=Xg^=0 

iflien f. |\e, 3 = 0 , 1 ,...., q, for some prescribed values of e, 
the process is considered to have converged and the parameter 
estimates are obtained from Eqn. 3.41. 

Initial Estimates of the Identified Models ; The initial 
estimates of the parameters were obtained by the methods 
discussed above. The forai of the identified models with initial 
estimates and the residual variance are shown in Table 3.2. 

Maximum Likelihood Estimation of Parameters : There are two 
methods for estimating the parameters of ARIMA models. They 
are 

(i) Estimation using Bayes' theorem 
( 11 ) Estimation using Maximum likelihood principle. 

In the present study, the latter method suggested by Box and 
Jenkins [65] is used for estimation of the parameters. For a 
sample of N observations z(t), there will be a probability 
distribution pCz/"^') depending upon some unknown parameters fe. 
The vector refers to p+q+1 parameters ( ©, o^) of the 

ARIMii model. Before the data are available, p(z/^^) will 
associate a density with each different outcome z(t), for 
fixed ^ . After the data have become available, there will be 

various values of ^ which might have given rise to the fixed 
set of observations z(t). The appropriate function for this 
purpose IS called the likelihood function 1 (^/ 2 ), which is 
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of the same form as p(s/’'^*), but in which z(t) is now fixed and 
~ IS variable. It is more convenient to work with the log- 
likelihood function l("^V2) = In . 

The unconditional log-likelihood function corresponding 
to N observations assumed to be generated by an ARIMA model 
IS given by 

1(9, 9,0 ) = f (9,9) - N In a - ( 3 . 44 ) 

e _ - ^ 2 

e 

where f(9,9) is a fimction of 9 and 0 , The unconditional 
sum of squares function 8(9,9) is given by 

K i 

3(9,9) = y [£(t) |9,9 ,z]^ ( 3 . 45 ) 

t=-co 

Usually f(9,9) is important only for small H. For moderate 
and large values of U, Eqn. 3.44 is dominated by S(9,9)/2 ct^ 
and thus contours of the unconditional sum of squares function 
in the space of the parameters ( 9» ® ) ^.re very nearly 

contours of log likelihood. Thus the least square estimates 
obtained by minimizing sum of squares given by Eqn. 3 . 45 , will 
provide very close approximations to the maximum likelihood 
estimates. 

Calculation of unconditional sum of squares; Consider the 
general ABMii model 

9(B) 5 (t) = 0 (B) e(t) ( 3 . 11 ) 

The above equation can also be written in terms of backward 
model as 
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9 (F) x(t) = &(F) e(t) (3.46) 

where P is the forward shift operator defined by 

F x(t) = X (t+l) and x(t) = x(t+k) (3.47) 

Taking conditional expectations of Bqns. 3.46 and 3.11, that is, 

9 (F)[x(t)] = 0(F)[e(t)] (3.48) 

9 (B) [£(t)] = e(B) [e(t)] (3.49) 


Equation 3.48 is first used to compute back forecasts and then 
Eqn. 3.49 is used to generate the [e(t)]*s. If the forecasts 
are found to be negligible in magnitude beyond some lead time ol, 
the recursive calculation goes forward with 

[e(-3) \ 9,0,2] =0 0 = 0,1,2,... 

[«(-3) I t,©:, 2] =0 3 > Q-1 (3.50) 

The uncondition sum of squares is then calculated using 

3(9,0) = [£(t)]2 (3.51) 

The least squares estimates of the parameters are obtained by 
minimizing che sum of squares given by Eqn. 3.51. For a 
purely autoregressive process residual e(t) is linear in the 
autoregressive parameters 9*3; however, for a purely moving 
average process e(t) is nonlinear function of parameters. The 
linearization of the model is done using Taylor series, let 
represents p+q parameters (9, 0). Expanding [e(t)] in a 
Taylor series about it^s value corresponding to the guessed set 
of parameter values ^2,o» '*■ ^(p+q) 0 i 
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gets 


and 


p+q 


E(t)=£(t) - N ('^ - ? )X 

Z__ XfO' x,‘ 

x^l 

where e^(t) = [£(t) j z, ] 


(3.52) 


o 


1 > t -A s • 

>» X 


X O' 


i ^ -»b 


§= -- 


(3.53) 


If X' IS the (N+Q)x (p+q) matrix of x , , the (N+Q) equations 

- i,x 

( 3 . 52 ) may be expressed as 


£.(t) = X' (5- ij + e(t) 


--0 


- O' 


(3.54) 


where eQ('t) a-nd e(t) are column vectors with (N+Q) elements. 

The adjustmentsC J - ^), which minimize S( 4 ) = 3(9,0) 

= (e_(t)) (e(t)) IS then obtained by linear least squares. The 

values of parameters which minimize residual sum of squares are 
obtained by the constrained optimization method suggested by 
Marquardt [71»'72]. The computer algorithm for this method [ 73 ] 
IS described in Appendix D. The residual variance is then 
estimated using 


ja _ S(t} 


e U->-q 

and the variance-covariance matrix by 


( 3 . 55 ) 


I'li) = i (5-56) 

The standard error of the parameters is obtained by taking 
square root of the diagonal elements of the variance-covariance 
matrix. The 95 per cent confidence limit for is given by 
+ 1.96 times the standard error of 
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Tile maximum likelihood estimates of the identified 
model were obtained by the methods discussed above. The summary 
of the models fitted to the series A to P is given in Table 3.3. 
Also shown in the t'able are the standard error of the parameters 
and the residual variance. 

The aiialysis of multiple tine series (Chapter 47/' 
requires that the sample size of all variables should be equal. 
Prom Table 3.3 it can be seen that a model (2,0,0) has been 
fitted to 270 observations of hot metal temperature. The miss- 
ing values were calculated using the equation 

y(t) = 9 ^y(t-l) + ^'‘^y(t-2) + (l-9^-92)^i 
and the paramecers of the model were reestimated using 350 
values. The process of calculating missing values and estimating 
the parameters was repeated until there was no change in the 
parameter estimates in two successive iterations. The final 
parameter estimates of the model for hot metal temperature 
are also given in Table 3.3. 

3.5.5 Validation of the Model ; 

The model having been identified and parameters esti- 
mated, diagnostic checks are then applied to residuals to check 
the adequacy of the fitted model. The fitted model is considered 
to be adequate if the residual series constitute a.n independent 
normally distributed series. Normality of residuals is tested 
by Chisquare rest and serial independence is tested by 
correlogram analysis and spectral analysis. 
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Formality of Residuals ; The samole space is divided into I 
mutually exclusive classes with a class frequency of 5 or more. 
Let p(i) he the proha.hility that the variable belongs to ith 
class and if e(i) and e(i+l) ere the limits of the ith class 
int erval then 

p(i) = F[e(i+1)] - P|[e(i)] (3.57) 


The value of P[e(i)] is calculated from the Formal probability 
distribution table, let f(i) be the observed frequency of the 
sample from the ith group. The Chi-square statistic < is given 


Fp(i) 


(3.58) 


If J IS the number of parameters estimated, then theoretically 

has a chisquare distribution with (l-J-l) degrees of freedom. 
Let '\^^(a') bo the value of a' per cent confidence level 

for the above degrees of freedom as obtained from statistical 
taoles. If the calculated value of is less than the theo- 

retical value, the residaals can be assumed to be normally 
distributed. 


Gorrelogram Analysis ; The autocorrelation function of a pure 
random process is uncorrelated and distributed approximately 
normally about zero with variance l/F and has a standard 
error of (l/F)^ . The autocorrelation function of the residual 
series is plotted. If a.11 autocorrelation coefficients lie 
within 95 per cent confidence limit (+ 1.96 times the standard 
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error) thi syndicates that the residual series is pure random 
and the model is considered to he adequate. However, it was 
pointed by IHirbin (see [65]) that using H ^ as standard error 
for residual autocorrelation function may underestimate the 
significance of apparent discrepencies. Hence, rather than 
considering the residual autocorrelation coefficients indivi- 
dually, it IS desirable to consider the smallness of first K 
autocorrelations. This is done by calculating a Q-statistic 
given by , 


Q = H 


1=1 


r^(e) 


(3.59) 


The fitted model is considered to be adequate if Q is less 
than the value of '|^^(a) at a* per cent confidence limit at 

i 

(K-p-q) degrees of freedom. 


Spectral Analysis ; The test for independence of residuals 
can &lso be performed in the frequency domain. The spectral 
density function for a pure random (white noise) process, is 
constant over the frequency range 0 < f f jand is equal to 
variance of the process, which is sa.me as the average of the 
computed values of spectra. The sample spectral density 
function IS distributed about its population value according 
to a distribution, where X- is the number of degrees 

of freedom given by 


L> = 


2N 

M 


2 

3 


( 3 . 60 ) 
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where M stands for the number of lags estimated autocorrelation 
coefficients are used for spectral estimation, From the 
Ghisquare Tables the values of ) and 3--^ ci:/2,i‘) 

are read at a given significance level a. From these values 
i are determined. If G-(f) denotes the spectral density 
for pure whitenoise and G^(f) that for given sample then 
Gj^(f)/G(f) has / JJ distribution and the series is not 
significantly different from pure random series at confidence 
level (l-a) if 


(1 - a/2. V) 
V 


4 


G, (f ) 


j^(a/2, V) 

2J 


(3.61) 


Validation of the Models Fitted to the Plant Data ; The 
autocorrelation functicms were calculated for residual series 
obtained from all models and the Q-statistic given by Eq_n. 3.59 
was calculated. The calceilated and theoretical valuesoof 
chisquare have been given in Table 3.4. The simplest model 
for which the calcijlated value is less than the theoretical 
value at 95 per cent confidence level, was selected for each 
series. 


The normality of the residuals was tested by chisquare 
test given by Eqn. 3.58. The calculated and theoretical 
value of chisquare at 95 per cent confidence level have been 
given in Table 3.5. It indicates that all residuals are 
nearly normally distributed. 



TABLE 3.4 


TEST FOR SERIAL iroSPEFDERCE OF THE UNIVARIATE RESIDUALS 


Series 

Model 

Q-Statistic 


Remarks 

No. of degrees Theoretical 
of freedom value at 95 

Calculated 

value 



per cent CL 



1 

2 

3 4 

5 



(1,0,0) 

11 

19.70 

25.30 

Model 


23 

35.20 

35.67 

rejected 


29 

42.00 

37.47 

(2,0,0) 

10 

18.30 

24.32 

Model 

22 

33.90 

34.87 

rejected 


28 

41.30 

36.64 


(3,0,0) 

9 

16.90 

14.41 

Model 


21 

32.70 

25.42 

accepted 


27 

40.10 

27.25 


(1,0,0) 

11 

19.70 

20.42 

Model 

23 

35.20 

36.66 

rejected 


29 

42.00 

41.54 


(2,0,0) 

10 

18.30 

18.52 

Model 


22 

33.90 

34.70 

rej ected 


28 

41.30 

39.57 

(3,0,0) 

9 

16.90 

11.60 

Model 


21 

32.70 

22.55 

accepted 

• 

27 

40.10 

27.79 


(1,0,0) 24 36.40 14.11 Model 

accepted 


(1,0,0) 

11 

23 

29 

19.70 

35.20 

42.60 

19*60 models 

on rV are adequate 

20-Si butd.oL) 

(0,0,1) 

11 

23 

29 

19.70 

35.20 

42.60 

g CO selected 

13*64 

iR*n? correlogram 
and partial- 

(1,0,1) 

10 

18.30 

5.01 correlogram 


22 

33.90 

10.49 indicate that 

(1,0,2) 

28 

41.30 

11.65 the model 

9 

16.90 

4 . 60 should have 


21 

32.70 

10.34 both AR and 


27 

40.10 

11.43 MA terms. 
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Table 3.4 (contd) 


“■ 1 

2 

3 

4 

5 

5 “ 


(3,0,0) 

9 

16.90 

15.48 

Model 

J 


21 

32.70 

21.31 

accepted 



27 

40.10 

22.95 


D 

(2,0,0) 

10 

18.30 

8.75 

Model 



22 

33.90 

21.11 

accepted 



28 

41.30 

27.48 


E 

(1,0,0) 

11 

19.70 

10.18 

Model 



25 

35.20 

20.42 

accepted 



29 

42.60 

23.97 



(1,0,0) 

11 

19.70 

26.06 

Model 



23 

35.20 

41.47 

rejected 



29 

42.60 

52.35 



(1,0,1) 

10 

18.30 

17.54 

Model 



22 

33.90 

32.65 

rejected 



28 

41.33 

42.20 


P 







(1,0,2) 

9 

16.90 

16. 65 

Model 



21 

32.70 

30.92 

rej ec ted 



27 

40.10 

41.45 



( 2 , 0 , 1 ) 


9 

21 

27 


16.90 

32.70 

40.10 


15.78 

27.67 

37.75 


Model 

accepted 
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TABLE 3.5 

TEST FOR NOR]>L'LITy OF UNIViiRlATE RESIDUALS 


Series 

Degrees of 
Freedom 

Theoretical 
Value at 95 
per cent CL 

Calculated 
Va lue , 

Eqn. (3.58) 

.,(t) 

18 

28.87 

27.72 

e2(t) 

25 

37.65 

29.06 


29 

42.56 

37.80 

e4(t) 

20 

31.41 

28.41 


25 

37.65 

36.78 

egCt) 

25 

37.65 

32.52 
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The residual autocorrelation functions for fitted 
models have been plotted in Fi^s. 3.19a to 3.19c. The power 
spectra for tne fitted models have been plotted in Figures 
3.20 to 3.25. Also shown in these figures are the 95 per cent 
confidence levels, calculated by Eqn. 3.26 for autocorrelation 
function and Eqn. 3.61 for spectral density function. In 
calculating the latter, the series were standardized to zero 
mean and unit standard deviation. Hence, the theoretical value 
of spectra which is equal to the variance of the series is 
1.0. All these results support the fact tnat the residuals 
series are not different from the pure white noise series at 
95 per cent confidence level. 

Another test for adequacy of the fitted model is to 
calculate power spectra of the original series using the 
estimated AR and MA parameters. For a general ARIIA (p,q) process, 
the power spectra is given by [74,75] 


0(f) = 2a^ 
e 



t 

j 





g-0 2q7tf 

a ^ 

-0 2p7tf I 

(3.62) 


2 


The above equation can be simplified for different values o 
p and q. Power spectra for some processes are as follows; 
Second order AR process; 

2o^ 

_e (3.63) 

[l+cp ^+92 “ 29^(1-92)^002711 - 292 Cos 47 if] 
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I 



Fig 3 24 - Power spectra of univanate residuals of silicon 

content senes E 



Frequency f m cycles per cast 
Fig 3 25 -Power spectra of univariate residuals of sulphur 
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Third order aR process: 

2a^ 

G(f) = — 

[ (l-9^Cos2'n:f - cp2^os4Tt:f - 9^003611:1)^ + 

(9^Sin2irf + 92Sin4ii:f + 9^Sin6ixf)^ ] 

(3.64) 


Mixed AEMA (1,0, l) Process: 

l +©2 _ 2 © Cos2itf 

G(f) = 2cr^ ^ (3.65) 

- 29 ^ Cos27i;f 

The confidence limit for power spectra is calculated 
using Eqn. 3.6l, in which now G(f) stands for power spectra 
at a given frequency f of the series under consideration. The 
model IS considered to be adequate at 95 pei* cent confidence 
level if power spectra generally lies within the confidence 
limit. 

Using the models of best fit, power spectra and 
confidence limits were calculated for all series. The power 
spectra calculated using Eqn, 3.62 and confidence limits are 
also shi^wn in Figures 3.13 to 3.18. In all the cases, the 
spectra of series lies within confidence limit which support 
the adequacy of the fitted models. The sample spectra is nea.rer 
the lower confidence limit than the upper one. This may be 
due to the maximum likelihood estimation procedure adopted in 
the studjr and other procedures for parameter estimation may 
lead to different results. In this study, the fitted models 
were considered to be adequate. 
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3.5.6 Stability of Univariate Models : 

The stability of the fitted univariate models was 
tested by finding the roots of the characteristic equation 
3.21. The characteristic equations and the corresponding 
roots have been given in Table 3.6. In all the cases, it 
was found that the roots of the characteristic equations lie 
within the unit circle. Hence ell the fitted models are 
stable. 



1D3 


TABLE 3.6 

STABILITY OP IffilVARIATE MODELS 


Senes 

Order of 
the 
model 

Characteristic Equation 

Roots 

A 

(3,0,0) 

z^-0.977z^+0.197z-0. 151=0 

0.940; 

0.0185+0. <3 

B 

(3,0,0) 

z ^-0 . 694z^+0 . 004Z-0 . 121=0 

0.855; 

-0.0805±0.037d 

n 

(1,0,0) 

z -0.897 = 0 

0.897 

C\J 

o 

(1,0,1) 

z -0.282 = 0 

0.282 

”^3 

(3,0,0) 

z^-0.684z^4-0.205z-0.344 = 0 

0.890; 

-0.103 +0.613 

D 

(2,0,0) 

z^-0.237z - 0.308 = 0 

0.687; -0.45 

E 

(1,0,0) 

z -0.473 = 0 

0.473 

P 

(2,0,1) 

z^+0.222z - 0.354 = 0 

0.716; -0.495 


4"r-l-4-r -j— h-H" 



CHAPTER 4 


MULTIVARIATE TIME SERIES ANALYSIS 

The sequence of values of one variable constitutes a 
univariate time series. The sequences of values of a number 
of variables constitute a multivariate time series. The 
modelling and analysis of multivariate time series should ta.ke 
into account not only the serial dependence within each of the 
series but also the mutual or cross-correlation among the time 
series. 

4.1 GENERAL MULTIVARIATE MODEL: 

4.1.1 Autoregressive Moving Average (ARMA) Model: 

Let the multivariate vector of the variables be given by 

x(t) = [x^(t), x^Ct), jXj^Ct)]'^ 

where K is the total number of variables and T stands for the 
transpose. To represent the relationship among the time 
series, a general (p,q) order multivariate model of the follow- 
ing form may be considered: 

x(t+l) = ^2 x(t) + S.C'*^+l-p) + 1]_ B(t+l) 

....... + B E (t+2-q) (4.1) 

=q- 

where A_^ , i = l,2,....,p are the autoregressive (AR) coeffi- 
cients matrices and , 3=1, 2,..., q are the moving average (MA) 
coefficient matrices, all of dimension (KxK) and E(t) stands 
for the vector of residual series. Equation 4.1 can be consi- 
dered as an extension of the univariate ARMA model to multivariate 
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domain. With. adeq.aate data, this equation can be solved for 
all unknowrns. Since there are large number of parameters, the 
result may be data dependent or unstable. 

4.1.2 Triangular Two Sided Moving Average (TTSMa) Model; 

Since en autoregressive model can be expressed as an 
infinite order moving average model, the general ARIIA model can 
also be expressed as an infinite order moving average model. 
Considering finite moving average terms, Phadke et al.[2] have 
proposed a triangialar two sided moving average (TTSWA) model, 
which may be represented as follows: 



+D _2_il(‘fc+l) + + * • • • • 

5n 1 (4.2) 

where p(t) is a vector of serially and mutually independent 
r e sidua 1 s eri e s and , i — —q^ , ,...*, 0 ,l,..o.,q 2 are lower 

triangular matrices. Equation 4.2 cen also be written in terms 
of backshift operator B as follows: 


x(t) = D ^ B“'^l r)(t) + D B“V^ p(t) + +D -,B"^ p(t) 


+D 

=0 


rj(t) + B r|(t) + 


.+D 

= 0 . 2 “^ 


B 


12-1 


Ti(b)+D B'lppCt) 

-l2 - 


(4.3) 

or x(t) =hpi'(B) rj_(t) (4.4) 

where (B) is a lower triangular matrix of the form 
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and 



T 

ri(t) = [ri^Ct), ri2(t) jrig-(t)] is a vector of K serially 

and mutually independent series. 


4.2. DECOUPLED MULTIYARIATE MODEL; 

In developing a decoupled multivariate model, initially 
a univariate stochastic model is fitted to each of the series, 
thus accounting for the internal correlation. The resulting 
residuals, which are serially independent and normally distri- 
buted are called ’prewhitened series! A multivariate model is 
then fitted to the univariate residual series in terms of 
serially and mutually uncorrelated random components to explain 
external correlation. This approach for modelling of multivariate 
time series is advantageous because of the following reasons; 

(i) Since serially independent random variates are used 
in the multivariate model, a better performance of the model 
can be expected than with the serially correlated data; 

(ii) Since the parameters are estimated in stages, the 
number of parameters in the two stages are comparable respect- 
ively to that of univariate and multivariate models and the 
problem of simultaneously estimating all parameters is avoided. 


and 
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(ill) This procedure enables the most appropriate form 
and order of the univariate model to be chosen independently 
for each variable. The form and order of multivariate model 
are independent of tne univariate models. Thus, it affords 
greater flexibility in the choice of appropriate but different 
models for the univariate and multivariate time series. 


4-2.1. ABMA Model; 

Recently Ramaseshan et al. [76] and Erishnasami [77] 
have proposed a decoupled multivariate model. According to this 
method, let_e_(t) be a vector of serially independent random 
variables, obtained as residuals from univariate tine series 
analysis. All residuals have mean zero and standard deviation 
^ , . The set of series e, (t) is then standardized to a E(t) 
senes with mean zero and unit standard deviation. Hence, 


S(t) 

E(t) = [-1 


A 

a 




a 


■1 


ek 


B(t) constitutes the multivariate series of serially independent 
components and is related to the serially and mutually uncorre- 
lated series by a multivariate ARMA model of order (p,<l) , as 
follows 

E(t) 




^ 1 (t-n) +'S p(t+l~v) (4.6) 


If u=v=l, Eqn. 4.6 becomes 

E(t) = E(t-l) + p(t) (4.7) 

let and M ^ be lag zero and lag one correlation matrices. 
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In order to preserve these correlations the coefficients of the 
model given hy Eqn. ^.7 can he estimated from 

=1 " i -1 

and £liI=So-Slii (4.9) 

This IS essentielly similar to the method of moments. The 
matrix is a lower triangular matrix. In this method all 
elements of matrix and all lower triangular elements of the 
matrix are to he estimated whether they are statistically 
significant or not. However, some of the elements may he 
statistically insignificant and it is not necessary to estimate 
all elements, in such case. Secondly Eqns. 4.8 and 4.9 are used 
when the multivariate modelis first order autoregressive and 
first order moving average. If the order of the model is higher 
than (l,l), the estimation of parameters becomes more complicated 
because a large number of them have to he estimated and moreover, 
these estimates may not he reliable. 

4.2.2 TTSriA Model ; 

As in case of general multivariate model, the multi- 
variate ARMA model given by Eqn. 4.6 may he expressed as an 
infinite order moving average model. Phadke et al. [2] have 
developed a TTSMA model of finite moving average order which 
IS given by 

e(t) = 2 (t) (4.10) 

where e(t) = [e^(t), ] (i-S the vector of 
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univariate residuals and is the lower triangular matrix 

of the form given hy Eqn. 4.5. The advantage of this method 
over previous method is that all elements of the lower tri- 
angular matrix need not have to he estimated. The 

statistical significance of the eletonts is obtained by plott- 
ing cross-correlation function betxreen e(t) and pCt). Hence, 
only those elements which are significant have to be estimated, 
say, by the method of least squares. 

The block diagram of Eqn. h.io is given in Figure 4.1. 

Thus, 

ei(t) = ■n^(t) (4,11) 

The series S 2 ('fc) consists of two parts, 62 3 _('t) ^2 

£2(t) = £2 1*^"^^ ^2 2^^^ (4.12) 

Out of these, £2 ^('^) the projection of e 2 ('b) series on the 
space of P 2 _(t) series and has the model 

The other component £2 2 ^’^^ orthogonal to series and 

has the model 


^2.2^'^^ = *422(3) Ti2('^^ 


(4.14) 


Thus, 


(t) =4-^^(B) p^(t) + 4 -^p 2(®^ Bp('t) (4.15) 


'1 


22 


In general, the series £ (t) consists of j orthogonal components, 
£^ k(^) ” 4^3 k^®^ Ei^('t)» =1,2,. ..,j. The first j-l compo- 

nents stand for the projection of £ (t) on the space of 

t] 
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T], (t), k=l,2 , . . . , 3 - 1 . The 3 th component e (t) defines 
a new whitenoise series ortnogonal to > 

k=l,2 3 -1. Hence, this-model ^represents Gram-Schmidt 
orthogonalization of stationary multivariate time series £(t) 
into K orthogonal vectors ) j k=l, 2 ,....K. 

The method adopted by Phadke et al to identify and estimate 
the parameters of the matrix ^(B) is as follows; 

( 1 ) Since residual series obtained by uni- 

variate modelling of the series, x^(t), it is a prewhitened 
series. Hence, the Eqn. 4.11 takes the form 

s^(t) = p^(t) (4.16) 


( 11 ) identified by ploxting the cross-correlation 

function between £ 2 ^'*^) Apnroximate estimates of 

are obtained from 


U 

in which r 
between £2 
calculated 
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(B) = (Yar eg(t)) = 


(Var p^(t) 


(4.17) 


' 2^1 


(B) stands for the cross-correlation function 


(t) and p^(t) series at lag B. £2 ^(t) is then 
from Eqn. 4.13. 


Ej.lCt) = %i(B) T,^(t) 


( 4 . 13 ) 


(ill) On subtracting £2 from e 2 (t) gives the series 

82 2 ^'^) according to Eqn. 4.12. 

(iv) £2 2 ^"^^ expressed as Eqb. 4.14 and the paramexers 
of "^ 22 ^®^ estimated by univariate modelling of time series. 
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This procedure is repeated for all series. That is, 

(a) Calculate cross-correlation function between e (t) and p (t), 
j = 1,2,..., 1 - 1 ; and hence identify 

(h) Multipljr £^(t) by corresponding 3=1,..., i-l to 

give the (i-l) components of £^(t) which are nothing but pro- 
jections of £^(t) onp^(-c), 3 = 1 , 2 , . . . , 1 - 1 . 

(c) Subtract these components from £^(t) to give the series 

£ (t) which IS orthogonal to p (t), 3 = 1 , 2 , ...i-l. 

-L _L J 

(d) Identify and estimates by univariate modelling 

of £ (t) series. 

11 

(e) Finally estimate the parameters of (B), 3=l,2,...i-l by 

^ 3 

the method of least squares. 

Tn the above method, the Grg,m-Schmidt orthogonalization 

orocess has not been used directly. Secondly, in identifying 

(B) , the cross-correlation function was calculated between 
13 

£^(t)^and p^(t), IloWjS^Ct) consists of i components. In 
particular, the presence of noise p^(t) or (B)p^ ( t ) 

corrupts the cross-correlation function between £ (t) and p (t). 
Hence, in the present study, the above method of developing the 
matrix 'M^(B) has been modified. 

4.2,3 Generalized Method for Parameter Estimation ; 

The method of building the multivariate model can also be 
generalized which leads to a simplified method for parameter 
estimation. From the set of vectors e(t) one can obtain a set 
of orthogonal vectors §.(t) by means of Gram-Schmidt orthogo- 
nalization. The algorithm [ 78] is as follows, and details of 
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which are given in Appendix B , 


p^(t) = e^(t) 


(4.18) 


a^(t ) 


p^(t) 


PT(t) 

||Pl(t)|| 


e (t) 

J 


(a (t)*£ (t)) a (t) 

I X J X 


1< 3 K 


P, (t) 
(t) = 

inxitT 


p^(t) 


11^1 


(4.19) 


( 4 . 20 ) 


( 4 . 21 ) 


in which ci:^(‘fc) represents ixh orthonormal vector and ljp^(‘t)|< 

stands for the nora of the vector p (t), (a (t)*e.(t)) stands 

i 13 

for the inner product of a. (t) and e (t). Fow writing Eqn. 

3- 3 

4.20 for , 


P,(t) = e„(t) - (a^ (t ) •c,(t) ) a,(t) 


( 4 . 22 ) 


or “ ‘^ 12 ‘“l^^^ 

where = (Q:^(t ) • £2 ( t ) ) 


(4.23) 


Substituting for a^(t) from Eqn. 4.19 ir. 4.23 

P2(t) = £2^^^ “ ^^12/-I1^1^'^MP ^1^^^ 

which may be written as 

e2(t) = H2^ P^(t) + 


(4.24) 


(4.25) 


where X 


(4.26) 


Similar Ijr, 
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p^(t) = e^(t) - (a^(t) • £^(t) ) 


a^(t) - (oc^Ct) *£^(1)) a^{±) 

(4.27) 


which when combined with Eqn. 4.21 becomes 


^^{t) - ej(t) - ~ iipffnn ^ 2*'*^' 

' 'i (4.28) 


where 

‘^13 “ (o^qC'fc) '^^(■t) ) 



623 = (a 2 (t ) -sj (t) ) 


or 

e^(t) = p^(t) + H ^2 ^ 2 ^^^ ^ 3 ^^^ 

(4.29) 

where 

H 31 = 6 ^^/ j|B^(t)|j 



Hj 2 = 623/l|P2(t)|| 

(4.30) 


In general, 


£ (■(:) 
J 




where H 


2'^ 


(a^(t) -e^ (t) )/ 


= 6 , 


k3 




(4.31) 


(4.32) 


The Eqn. 4.31 can be written for all va.riables from o=l,2,...K 
as follows; 


\(t) ' 

i 

1 




Pi(t) 

£2 

i 

H21 

1 



PqCt) 

£^(t) 

: 

= ' 

H . 

3I 

H^. 2 ... 



P^Ct) 

ieu{t) 

} i 

R 

H 

H_ .. 1 


Pg(t) 1 

1 K1 

K2 * " 

ii*n-± j 

L 


(4.33) 
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or ^(t) = H p(t) (4.34) 

This IS known as a principal component model. 

Since g.(t) IS a set of orthogonal vectors, the dot 

product between any two series p (t) and p (t) , 17^3 will be 

zero. Hence, cross-correlation function at lag zero will be 

zero, however, cross-correlation function between P^(t) and 

P (t) may be significant at lags other than zero. Let r„ „ (k) 

3 P2P1 

be the cross-correlation coefficient between series P2('*^) 

P^(t) series at lag k. From the cross-correlogram, the 

coefficients which s-re significant at 95 per cent confidence 

level are determined. Let r^ ^ (B) represent the significant 

cross-correlation coefficients at lag B, then a quantity 

IS calculated from 

= T ( 4 . 35 ) 

(Var p^(t))2 

e IS then calculated from 

^(t) = H 2 ^-p^(t) + (4.36) 

and IS subtracted from give 2 ^'^^ series, 

IS then expressed as Eqn. 4.14 and the parameters of ^22^®^ 
are estimated by uniAT'ariate modelling of time series. Thus 
e2(t) can be expressed as 

e2(t) = H2^P^(t) + A23_(B) p^(t) + U^22^®^ ^2^^^ 

= (H2^+X2 i(B))Pi(^) + ^2^^^ (4.37) 
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The parameters of the model given by Eqn. 4.37 are estimated 
by multinle linear regression analysis. The method of centra- 
lization and correlation matrix [79] has been used for 
regression analysis. This procedure is repeated successively 
for all series. 


4.3. TRAIT SEER OTGTIOR MODEL; 

The decoupled multivariate model developed by Phadke et 
al. IS useful in determing the transfer function model for a 
physical system under control. Consider a stable, linear, 
closed loop system which has m-manipulatable variables, X (t), 
0 = 1 ,..., m, and n-output variables Y^(t), i=l,....n. Then 
letting X(t) and Y(t) to represent corresponding m-dimensional 
input vector and n-dimensionaO output vector, the model for 
the plant can be written as 


Y(t) = Y(B) X(t) + N(t) 


(4.38) 


where V(B) is a nxm matrix wnich represents plant dynamics. 

Every element 7 (B) of V(B) represents transfer function 

10 = 

between the output Y. (t) and the input X (t), and has the 

^ 0 

form 


10 w (0) + w (1) E+....+ M ( 3 ^ ) B®io 




B 


-IJ 


11 


10 11 . 


{4.38n 


1] ij 

where b^^ represents the time delay or dead time between the 


innut output Y^(t). 
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N(t) represents the plant disturbance, which is n-dimensional , 
nonsingular stationary process and can be described by 

N(t) = S (■t) (4.39) 

where (t) is n-dinensional white noise process such that 
its variance-covanance matrix S[ ( (t))( v (t)) ] = 

IS a diagonal matrix. 

Let the characteristics of the feed back dynamics matrix C(B) , 
a mxn matrix, be similar to those of V(B); then the model for 
feedback can be written as 

X(t) = C(B) Y(t) + M(t) (4.40) 

where M(t) represents m-dimensional, nonsingular, stationary 
feed back disturbance and can be described by 

M(t) = *4^'(B) _V^t) (4.41) 

in which (t) represents m-dimensional white noise process 

such that its variance-covariance matrix is a diagonal matrix. 

If I and I represent mtn order end nth order unit 
n n 

matrices, respectively, then Eqns. 4.38 and 4.40 can be written 
as 


[I 

^=m 

- C(B) 

V(B)] 

X(t) = 

= M(t) + C(B) 

li(t) 

(4. 

.42) 

[I 

- V(B) 

C(B)] 

Y(t) = 

= V(B) M(t) + 

N(t) 

(4. 

.43) 

'-=n 




— — 




If 

£PB) 

= Cim 

- C(B) 





and 

D„(B) 

= [I 

- 7(B) 

o(B)r^ 


(4, 

.44) 


=2 

‘•=n 


— 




Then Eqns. 

4.42 £ 

'nd 4.4! 

5 become 
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X(t) 

= i^{B) K%) 

+ 5i(B) 

C(B) W\b) 

r'<t) 







(4.45) 

I(t) 

= S 2 

(B) V(B)y'(B) 

C*'(t) + 

D 

2 (B) 4^(b) 

lit) 







(4.46) 

Let 


II 






ill'®) 

= D^(B) "^''(B) 

< 

9 

(mxm) 

(4.47) 


Ii2<B) 

= D^(B) C(B) 

V(B) 

f 

(mxn) 

(4.48) 


isi'®) 

= D^CB) V(B) 

fw 

f 

(nxm) 

(4.49) 

and 

122(B) 

= 12(B) V(B) 

f 

(nxn) 

(4.50) 

Then ! 

Eqns. 4 

.45 and 4.46 can be written 

. as 



X(t) : 

= i^(B) r (t) 

^ ii 2 

(B) 

>'(t) 

(4.51) 


T(t) : 


^22 

(B) 

|'(t) 

(4.52) 

which 

may he 

combined to give 





Z(t) : 

= 1(B) u(t) 




(4.53) 

where 

z(t) 

= (x(t), i(t))^ 





and 

U('t) 

II 




The variance 

-covariance matrix, 

of 

white noise 

series 


^ It) IS the mxm principal minor of the variance- 

variance matrix of u(t), standing in the upper left corner. 

The nxn principal minor of standing in the lower right ■ 

corner is equal xo , the variance -covariance matrix of ' ^ (t' 
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The matrix 1(B) contains terms involving only positive povrers 
of B. However, the elements 3 > k of the matrix ^'(B) 

of TTSMA model may have poles outside as well as inside the unit 
circle. The matrix L(B) may be obtained from TTSWA model by 
removing poles and zeros inside the unit circle by using shift 
transformation given in [80 ]. If raw (original) data are used 
in the multivariate analysis, then Z(t) of Eq_n. 4.53 is equal 
to x(t) of general TTS'IA model given by Eqn, 4.3. If the 
prewhitened series are used for multivariate analysis, then 
Z(t) is equal to £(t) of Eqn. 4.10, 

The matrices Y(B) , C(B) , '-r (B) and i (B) are obtained 
from the matrix l(B) as follows; 

(i) Partition the matrix L(B) as follows; 


1(B) = 

hips) 

1 

? 

■ 


— — 

« « 

— 


L2i(B) 

1 

< 

122(B) 

-J 


where the dimensions of submatrices are given in Eqns. 4.47 
to 4.50. 

(ii) Combining Eqns. 4.49 and 4.50 

V(B) ='f^(B) LJ^(B) B23_(B) Y“^(B) (4.54) 

similarly Eqns. 4.47 and 4.48 give 

C(B) = ^f/’(B) ^}{B) (4.55) 

(ill) From Eqns. 4.50 and 4.44 one gets 
^4^^(B) = ^22^®^ - (4.56) 
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which on combining with Eqns. 4.54 and 4.55 gives 

^^(B) = ^2^^ " ill il2^®^ 

Similarly Eqns. 4.47, 4.44, 4.54 and 4.55 give 

= Lii(B) - B^2^B) L~^ (B) ^21^®) (4.58) 

Thus, from the multivariate model, various transfer functions 
(i.e., for plant, feed back, plant noise and feed backn noise), 
canbe obtained. 

4.4 ANALYTIC Al PROCEDURES; 

In order to identify and estimate the parameters of the 
multivariate model (Eqn. (4.10) it is necessary to calculate 
cross-correlation functions between the series e(t) and T).(t). 

Just as autocorrelation function and tne spectral density 
function are helpful to identify univariate time series models, 
the cross-correlation function and cross spectral density 
function are helpful in identification of the multivariate 
models. 

4.4.1 Cross-Correlation Function ; 

It determines the linear dependence between two series. 

The cross-covariance function between two series x(t) and y{t) 
at any lag k is given by 

Yj^yCk) = cov(x(t) 'yCt+k) ) = E[ (x(t )-p,^) (y (t+k) ) ] 

The theoretical cross correlation coefficient at lag k is defined 

By 



-L^± 


^ (]^) ^ ^ QJ (4.59) 

(var x(t)*var y(t+k))^ '^x '^y 


k=0 f +1 f +2 ...... 

Unlike the autocorrelation function, in general, cross- 
correlation function IS not symmetrical about k=o that is 
V (k) IS not necessarily equal to K (-k) . 


The cross-correletion function is usually estimated from 
N-k ^ M N-k 

H-k U"' - Tyrra > y 


-U X 


■T T 

(var x(i))2 (var y(n-k))^ 


where varx(i) = 


N-k 


N-k 

y x(i)^ 


(4.60) 

T N-k 

7 

(N-k)2 


N-k 
_ , 

1=1 


U-k 


var y(i+k) - ~ 


(iT-k) 


z _ 

1=1 


standard Error of Oross-Correlation Coefficients A crude check 

as to vrhether the cross-correlation coefficient ^ (k) is not 

xy 

different from zero is made by comparing the corresponding 
estimate of the cross-correlation coefficient with their approxi- 
mate standard errors obtained by a formula given by Bartlett. 

SE[r^y(k)] :Ai (N-k)"^ (4.61) 

which is used to test the significance of cross correlation 

cross- 

betxfeen two white noise processes. All ^correlation coefficients 
beyond + 1.96 times the standard error are considered to be 
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sigr.icant at 95 per cent confidence level. The two series 
are said to be mutually independent if the cross-correlation 
coefficients lie generally within + 1.96 times the standard 
error. 

4.4.2 Oross Spectrum ; 

The correlation between any two series in the frequency 
domain is represented by the cross spectrum. The sample cross 
spectrum contains two different types of information about the 
dependence between the two processes. It is a complex quantity 
and may be written as product of a real function called the 
sample cross amplitude spectrum and a complex function called 
sample phase spectrum. The cross amplitude spectrum shows 
whether frequency components of one series are associated with 
large or small amplitudes of the other series of the same 
frequency. The cross phase spectrum show's xfhether frequency 
components of one series lag or lead the components of the 
same frequency of the other series. The cross spectrum is 
obtained by taking Fourier transformation of the cross-covariance 
function. In the present study, only cross-correlation function 
IS used for identification purposes. 
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4.5 ANALYSIS AND DISCUSSION OP RESULTS; 

As series C is nonstationary (Sabsection 3.5.2) it was 
divided, into three parts and, after standardization, different 
models were fitted to the parts separately. The residuals 
oLtained from each part had zero mean and different variances. 
Hence, residuals of series C as a whole is nonstationary. In 
developing a multivariate model two ways may be adopted, 

(i) All the series may be divided into three parts as series C 
and fit a multivariate model to each of the parts, and (ii) 
Standardise residuals of each part of the series C so that 
series C as a whole hag zero mean and unit standard deviation. 
In the present study, the latter procedure is used. The 
standardized series C, along with residuals A,B,D,E and F 
constitute the multivariate time series. In multivariate time 
series analysis, the series were used in the order A, B,C,D,E 
and F. 

4.5.1 Generalized Method ; 

This method is described in Subsection 4. 2. 3. Table 4.1 
shows norms of the orthogonal vectors t ) , k=l,....,6; inner 
products of a^j-Ct) and ej('b); and the parameters of the 

principal component model (Eqn. 4.32). The autocorrelation 
functions of all the orthogonal vectors have been plotted in 
Figures 4.2(a) and 4.2(b). It is found that all orthogonal 
vectors are pure random series. In order to identify ,v\^^(B), 
1=1,..., 6 , 0=1,..., 6 , i> 0 » cross-correlation functions 
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TABLE 4.1 




PAEAJIETER ESTIMx'.TES OF PRINCIPAL 

MODEL 

o 

o 

o 

k 

3 

Norms of 

Q. (t) Inner product 

. =Vt)-eAt) 

1 

2 

7.06 

0.1786 

0.0252 

1 

3 

7.06 

3.9006 

0.5512 

2 

3 

12.00 

-0.9576 

-0.0798 

1 

4 

7.06 

0.9843 

0.1390 

2 

4 

12.00 

0.7777 

0.0648 

3 

4 

18.25 

-0.5932 

-0.0324 

1 

5 

7.06 

0.2548 

0.0360 

2 

5 

12.00 

2.0611 

0.1718 

3 

5 

18.25 

-5.5208 

-0.3021 

4 

5 

16.90 

0.4034 

0.0238 

1 

6 

7.06 

0.6053 

0.0855 

2 

6 

12.00 

-1.1507 

-0.0960 

3 

6 

18.25 

0.9659 

0.0528 

4 

6 

16.90 

-2.0189 

-0.1190 

5 

6 

14.96 

-7.5162 

-0.5000 
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"between all series of the vector §,(t) were calcula^ted. The 
cross-correlation functions have been plotted in Figures 4.3 
to 4.7. These plots indicate that the series P^(t), 1=1,,.. ,6 
are not mutually correlated. Since ^(t) is a vector of 
serially and mutually uncorrelated random series, it represents 
the multivariate residual vector, i.e., 


g.(t) = (4.62) 

The multivariate model can be written as. 


£]_(■*:) 


1.0 



- 

e2(‘t) 


0.0252 

1.0 



^^(■fc ) 


0.5512 

-0.0798 

1.0 




0.1390 

0.0648 

-0.0324 

1.0 

£5(t) 


0.0360 

0.1718 

-0.3021 

0.0238 1.0 

_eg(t) 1 


0.0855 

-0.0960 

0.0528 

-0.1190 -0.5 1.0 

1 


[Ti^(t), , ^^(t), Tl^(t), Pg(t)]^ 

(4.63) 


With a ^ 

= 0.143, 

6^ 2 3 , 

0.410, 

2 

= 0.953, 

^1 


ri2 


^3 



= 0.815, 

0 ^ = 

0.676, 


= 0.641 

h4 


^5 


^6 



The multivariate model given by Sqn. 4.63 is also a principal 
component model. 

4.5.2 Testing of Multivariate Residuals: 

The multivariate residuals should be distributed normally 
and they should be serially as well as mutually independent. 
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Fig. 4.3 -Cross correlogram betweerr fi'j (t) and succeeding orthogonai vectors. 
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formality ; The normality of the residuals was tested using 
Bqn. 5.58. The results are shown in Table 4.2, In all the 
cases the calculated value of chi -square statistic is less 
than the theoretical value and so it is concluded that all are 
nearly normally distributed. 

Independence Within the Time Series : The serial independence 
of multivariate residuals is tested from correlogram of the 
residuals. It was found earlier that g.(t) = ii{t) , (Eqn. 4.62). 

The correlogram of Pj^(t), i=l,2,...,6 is shown in Pigures 4.2(a) 
and 4.2(b). They indicate that all autocorrelation functions 
lie within the 95 per cent confidence limit. It is hence 
concluded that the multivariate residuals are not significantly 
different from pure random series at 95 per cent confidence 
level. 

Independence Among the Time Series: The mutual independence 

of the multivariate residuals is tested by cross-correlogram 
between each pair of residuals. The cross-correlation functions 
between series P^(t) and succeeding orthogonal vectors p^(t); 
i = 1,...,5; 3 = 1,...?6 have been shown in Figures 4.5 to 4.7. 
Also shown in the figures are the 95 per cent confidence limits. 
From these figures it is concluded that the orthogonal vectors 
and hence the multivariate residuals(as g.(t) = r].(t)) are 
mutually uncorrelated. 

A Q-statistic similar to that given by Eqn. 5.59 but in 
which the autocorrelation function is replaced by cross-correlatio 
function is used to test the mutual independence of the series, 
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TABLE 4-2 


TEST FOR NORMALITY OP I«JLTIVARIATE RESIDUALS 


Series 

Decrees of 
Freedom 

Theoretical 
Value at 95 
per cent CL 

Calculated 

Value, 

Eqn. (3.58) 


18 

28.87 

27.72 

n2(t) 

22 

33.93 

28.76 


28 

41.34 

35.65 


21 

32.67 

30.36 


24 

36.15 

25.98 

ugCt) 

25 

37.65 

21.22 
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viz. , 

K» 

Q = ^ K (4.64) 

i=-K’ ^ 3 

The calculated value of Q is compared with the theoretical 
value at 95 per cent confidence level. If the former is less 
than latter, the series under consideration can be considered 
to be mutually independent. 

The chi-square test for serial and mutual independence 
of the multivariate residuals are given in Tables 4.3 and 4.4 
respectively, and they show that the multivariate residuals 
are serially and mutually uncorrelated, 

4.5.3 Transfer Function Model ; 

The details of developing transfer function model have 
been given in Section 4.3. As TTSIM, model is also a principal 
component model, the vector u(t) of Eqn. 4.53 is equal to 
multivariate residuals vector T].(t) and the matrix 1(B) of 
Eqn. 4.53 is equal to the matrix ^t'■'(B) of. TTStl4. model. 

Thus, 


1.0 


i 




0.0252 

1.0 

1 




0.5512 

-0.0798 

1.0 ! 




0.1391 

0.0648 

-0.0324 i 

1.0 



0.0360 

0.1718 

-0.3021! 

0.0238 

1.0 


i 0.0855 

-0.0960 

0.0528 ' 

-0.1190 

-0.5000 

1.0 


(4.65) 
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TABLE 4.3 

TEST EOR SERIAL INDEPENDENCE OP THE MULTIVARIATE 

RESIDUALS 


Series 

Q - Statistic , Eqn. ( 3,59 ) 

^i(t) 

15.33 

29.59 

32.37 


12.43 

24.87 

31.51 


5.87 

17.75 

19.44 

N^Ct) 

8.84 

23.09 

26.83 

R5('t) 

6.63 

19.64 

23.94 

Ug(t) 

20.04 

34.75 

39.80 

Degrees of 
Preedom 

12 

24 

30 


Theoretical Value 
of at 95 per 
cent Cl 


21.03 


36,15 


43.77 
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TABLE 4.4 

TEST EOR I-riTUAL INDEPENDENCE OF OTLTIVARIATE RESIDUALS 


i 

2 


Q-Statiptic 

K * 

= 350 / r^ (k) 

k=-K’ 


1 

2 

34.01 

43.51 

51.18 

1 

3 

16.73 

35.89 

43.87 

1 

4 

18.67 

51.91 

65.10 

1 

5 

31.15 “ 

49.24 

65.44 

1 

6 

18.85 

31.33 

42.50 

2 

3 

33.81 

64. 64 

74.54 

2 

4 

14.20 

57.52 

63.32 

2 

5 

16.94 

45.99 

50.67 

2 

6 

23.85 

55.33 

77.81 

3 

4 

19.94 

48, 94 

62.32 

3 

5 

35.74 

51.57 

59.99 

3 

6 

15.46 

43.12 

61.78 

4 

5 

30.01 

58.32 

78.73 

4 

6 

35.99 

54.92 

67.37 

5 

6 

13.99 

45.89 

70.31 

Degrees of 
freedom 


24 

48 

60 

Theoretical Value 
of X ^ at 95 per 
cent CL 

36.15 

65.40 

79.08 
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and the variance-covariance matrix of uCt) is 
0.143 

0.410 

0.953 0 I (4.66) 

0.815 

0 0.676 

0.641 ! 

_j 

The matrix l(B) is then partitioned as indicated in Bqn. 4.65, 


thus 

! 1.0 


L^^(B) = I 0.0252 

1.0 


0.5512 

-0.0798 

1.0 

1 0 

0 

0 

= 0 

0 

0 

0 

0 

0 

j' 0.1391 

0.0648 

-0.0324 

I2:^B) = 1 0.0360 

0.1718 

-0.3021 

1 0.0855 

-0.0960 

0.0525 

ri.o 

i 0.0238 

1.0 


1-0.1190 

-0.5 

1.0 


The variance-covariance matrix, X.,,., of fhe white noise process 

S. 

5 (t) is given by 
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'X. . = 


0.615 


0 


0.676 0 

0.641 


( 4 . 68 ) 


The model for blast furnace noise is calculated using Eqns. 
4.39 and 4.57; hence. 


(4,69) 


The variance-covarxance matrix > ,,, of the white noise 


■%(t) ■ 


1.0 

^1 

ixM ■ 

HjCt) 

= 

0.0238 

1.0 1 

. 


1 


-0.1190 

-0.50 • . . 1.0 


t-U) 


process 'J(t) is given b3r 
0.143 


v"- 


0.410 


0 


0 


0.953 


(4.70) 


-4 


The model for feedback noise is calculated using Eqns. 4.41 
and 4.58; hence 


‘M^(t) ' 


M2(t) 

= 

M5(t) 



1.0 

0.0252 


1.0 



It) '■ 


>. li 

t^t) 

y 

ip) 


(4.71) 


The model for blast furnace dynamics is calculated using 
Eqn. 4.54; hence 

”0.1554 0.0622 -0.0324 

Y(B) = 0.1987 0,1476 -0.3021 [ (4.72 

0.0587 -0.0918 0.0528 
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The feedback dynamics is calculated using Eqn, (4. 55); hence 

0 0 0 

C(B) =000 (4.73) 

_ 0 0 0 

This indicates that controls of the order of one cast interval 
and above are not present in the system. 

The above relationships (from Eq;as. >4*65 to 4.73) are 
in terms of prewhitened input and output variables. In order 
to obtain the transfer function model in terms of actual input 
and output variables, the multivariate model given by Eqn. 4.63 
has to be first recoupled. Recoupling consists of substituting 
for the vector, e(t) of univariate residuals in terms of the 
vector y;(t) of the actual input-output variables using the 
univariate models. 

A siimmary of all univariate models fitted to series A to 
E has been given in Table 3.3 and the models of best fit are 
selected in Table 3.4. Referring these Tables, we may write 
expressions for e^(t), i=l,2,...,6 as follows: 

(lr0T-0.'97B + 0.197B2 - 0.151B^) y^(t) = e^(t) 

(1.0-0.692B + 0.004B2- 0.121B^) y2(t) = e2(t) 

-(1.0-0.684B + 0.205B2- 0.344B^) y^(t) = £^(t) (4.74) 

(1.0%0.237B - 0.308B2) y^(t) = e^^(t) 

(1.0-0. 473B) y^(t) = £^(t) 

(1.0+0.222B - 0.354 B^)yg(t) = (l.O+O. 582B)£g(t) 

where y. (t), i=l,...,6 stand for the standardized series in the 
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order A to F, and in the above equation, the (3,0,0) model 
for the last 127 observations of blast flow rate series has 
been used. Equations 4.74 may also be written in matrix notation 
as follows: 


yi(t).i 

1 

"(1.0-0.978+0.1978^ 

-Q.151B^)"^ 

y-('fc) ! 

(1.0-0. 692B+0.004B^ jO 



-0.121B^)~^ 

y^(^) 


(1.0-0. 684B+0.205B^ 



-0.344B^)-^ 

1 



(1.0-0. 23 7B-0 . 308B^ )“^ 

0 ^ 



(1.0-0.4738)“-^ 

|y.(t)| 

. (1.0-0. 582B) . 

L. o 1 I 

(i.0+0.222B-O.354B^) 


[e^ (t) ,e2('t) (t) ,£g(t) (4.75) 

or, y{t) = U(B) £(t) (4.76) 

On combining Eqn. 4.76 with Eqn. 4.63 the relationship between 
actual input-output variables and the multivariate residuals is 
obtained as follows; 

y(t) = U(B) ^ (B) rL(t) (4.77) 

or ^(t) = 1> (B) Ti(t) (4.78) 

Equation 4.78 is similar to Eqn. 4.53 and by the method described 
in Section 4.3 the corresponding transfer functions can be 
obtained. The model for blast furnace noise for actual input- 
output variables is then as follows; 
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Ni(t) 1 


1.0 


i 

1 

(1.0-0.237B-0.508B2) 


N2(t) 


0.0238 1.0 



Ti70-O.473B) (1.0-0.473B) 


N^Ct) 


-0. 1190-0. 0692B -0.50-0.291B 

1.0+0. 582B 

1 ! (1.0+0. 222B-0.354B2) ^^;_°^®*|P® 

T1.0+0.222B 

-0.354B2) 



(4.79) 


The variance -covariance matrix of is given by Sqn. 4.68. 

The model for feedback -noise for actual input-output variables 
is as followss 


M^(t) ' 


(1.0-0. 97B+0.197B^ Q 

-0.151 B^)"^ — 

M2(t) 


0.0252 1.0 


(1.0-0. 692B+0.004B^ (l. 0-0 . 692B+0 .004B‘^ | 

-0.121B^) . -0.1213^) 

M^t) 

\ 

\ 

! i 

1 

1 

0.5512 -0.0798 1.0 

^ (1.0-0. 684B+ T1.0-0./684B+ (1.0-0.684B + 

0.205B^-0.344B^) (0. 205B2-0.344B^ ) 0. 205B^-0.344B^)| 


X [C^(t), ^^(t)]^ (4.80) 

i4. 

The variance covariance matrix of ^(t) is given by Eqn. 4.70. 
The model for feedback dynamics is given by Eqn. 4.75 as before 
indicating absence of feed back control of order of one cast 
interval or larger. 
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The elemerits V. .(B), i = 1 to 3 and j=l to 3, of the 
blast furnace dynamics matrix 7(6) are as follows: 


7ii(B) = 7 „j3(B) = 

7 i 2 (B) = 7^33(6) = 

7^^(B) = = 

721(B) = 7 ^g(B) = 

^22^®) = = 
723(B) = VQg(B) = 

7^^(B) = \p(B) = 

7^2(2) = = 

^ 33 ( 6 ) = 7qj,(B) = 


0.1554 - 0.1509B + 0.02963B^-0.023lB^ 

1.0-Q.237B -0.308B2 

0.0622 - 0.0431B + 0.00026B^-0.00755B ^ 

1.0 - 0.237B - 0.3086^ 

-0. 0324+0. 0222B - 0.0066B^+0.0111B^ 

1.0 - 0.237B - 0.3086^ 

0.1S87-0.1937B + 0.0391B^-0.03B^ 

1.0 - 0.473B 

0.1477 + 0.01023B + 0.0006B^-0.0179B^ 

1.0-0.473 B 

- 0.3021 + 0.206B - 0.062B^ + 0.1040 B^ 

1.0 - 0.473 B 

0.05872 - 0.02162B-0.02150B^-0.00223B^-Q.0052B^ 

1.0 + 0.222B - 0.354 B^ 

-0. 0918+0. 00993B + 0.0366B^+0. 0108B^+0.0065B'^ 
1.0+0.222B - 0.354B^ 

0. 0528-0. 0045B - O.OlOlB^-0. 011 9B^ -0.01066"^ 
1.0+0.222B - 0.354B^ 


In the above expressions 7^^^ (B) represents the transfer 
function between sinter-to-coke ratio (Series A) and hot metal 
temperature (Series D) and so on. On simplifying the above 
expressions we get 
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- 0.075B - 0.1545 +(o! 45 ^ 1 . 0 ) * (o?687B-l.oT 


VbdCb) 






0.0116 


0.0234B - 0.0186 + (o?45^1.0) (0.687B-i7^ 


n 4 - n 0 ♦ 0 99 , 0 » 00137 

-0.036B + 0.0, -56 - (0.45B+I.O) (0.687B-1.0) 

0.06358^ + O.O516B + 0.518 + I' l ' ^ 


Tjj(B) = 0.0578B^ + 0.0787B 0.144 - x g:^^^.i.o) 




V^pCB) 




= -0. 228^-0. 334B-1. 14 


0.838 


(0.473B-1.0) 


0.0147B^ + 0.01545B + 0.112 + (g:g^§Li :o7 * ' lojillli- ' 

1 . 0 ) 

-0.01825B^-0.0421B - 0.182 -(o.*495B-1.0) l0*.%6B+1.0) 


“ 0.03B^ + 0.0525B - 0.1480 + (o?495B-1.0) ‘^(0.716^1.0) 


0.018 


The transfer function between the input and output given 

by Eq.n. 4.38a aasumes a pure delay of order b. . and a time 

. “ 3 

constant represented by a number of linear difference operators. 
The above equations indicate that the transfer function can be 
considered as weighted sum of a number of dela.y operators and 
time constants. It is known that all discrete models with 
finite inputs are stable and a comparison with differential 
system indicates that the forward and backward approximations 
to differential system should be so chosen that they are stable. 
Hence, we may consider the above representation to consist of 
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terms involving shift operator B which represent the time 
delay and the other terras which represent the time constant. 
For example, the transfer function between the input A and 
output F consists of the following; 

(i) a dead time of 2 cast intervals with a weightage 
of 0.0147 

(ii) a dead time of 1 C8.st interval with a weightage of 
0.01545 

(iii) a weightage of 0.112 with no delay 

(iv) a weightage of 0.0727 with a time constant T of 

°1 

0.495 cast interval, and 

(v) a weightage of 0.0116 with a time constant T of 

°2 

0.716 cast interval. 

The total weightage ¥, on the input is hence 0.2338 and the 
relative xfeights for the operators to are, respectively 
0.063, 0.0661, 0.479, 0.311, and '©•.0814. ' .HeMe, the time 
constant and time delay for the system are given by 


T„ = T„ VL + T 
c C-, 4 Co 5 


Tj^ = 2V7^ + ¥2 


A block, diagram representation of the input-ou.tput relation- 
ship is given in Figure 4.8. The total weight on the input 
and the fractional weights as well as the time constants and 
dead times have been calculated for all input-output pairs and 
are shown in Table 4.5. 
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Some workers have derived the transfer functions for 
the input-outmt pairs, blast temperature-silicon content and 
blast humidity-silicon content. These have been shown in 
Table 4.6. It can be seen from this table that, in general, 
a time constant of about 6 hours exists between blast 
temperature and silicon content. The response of silicon 
content to blast hmidity is very fast, the time constant is 
of the order of two hours, which is less than a cast interval. 

In the present study, blast temperature was not taken into 
account because when data were collected, it was kept at a 
constant level of lOOO^C, Castore et al. [62] have shown 
that a dead time. of 0.75 times cast interval and a time constant 
of 0.75 times cast interval exist between ore-to-coke ratio 
and silicon content of hot metal. 

For the present study, the data on output variables 
were available only at cast inteivals but the data on input 
variables were available at one hour tine inter\^al. Hence, 
the input va.riables were averaged between successive casts as 
a result of which a time delay of the order of half cast 
interval already existed between inputs and outputs. 

It can be seen from Table 4.5 that the total lag between 
the average input in a cast interval and output at the end of 
cast is of the order of 0.35 to 0.59 casts. Hence, the actual 
lag between input and output is of the order of 0.85 to 1.09 
casts. This is comparable to the results obtained by other 


workers. 
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Worker 


Blast humidity- 
silicon content 
of hot metal 


Blast temperature- 
silicon content of 
hot metal 


Wood[40] 


-0.0580 

‘l+2.8s 


l+6s 


Barhieri 
[see 40] 


-0.055 
1+3. 4s 


0.0025 
1+6. 2s 


Staib et al, 
[see 40] 


-0.023 

l+2s 


0.002 

1+14 s 


Gastore et al. 
[62] 


•- 0.8 

1+1 i 6s 


0^21 

l+3s 


Vidal et al._Q^Q2i , 0.007 

[ 58 ] + rx+Y:5^' 


0.004 

l+6s 


Present 

study 


0.0378B^ + 0.0787B + 0.144+ 

-0.0037 

(0.473B-1) 
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Results from Table 4.6 indicate the time constant of 
the order of 1.6 to 3.5 hours for blast humidity - silicon 
content transfer function. Since the cast interval is of the 
order of 2.5 to 3.0 hours, 1.09 cast intervals correspond to 
around 2.7 to 3.3 hours and this is in agreement with the 
earlier results. Vidal’s result indicates a time delay of 
8 hours on 25 per cent of the input corresponding to an average ^ 
time delay of about 2 hours. This may be compared with the 
results of this study, namely, 14.3 per cent weight on two 
time delays and 29.8 per cent weight on one time delay corres- 
ponding to an overall time delay of 0,584 cast interval or 1.5 
to 1.8 hours. Since the multivariate model gives greater 
details of the transformation and gives comparable results to 
one of the input-output pairs, the other results may be consi- 
dered to be reasonable in the absence of information to the 
contrary. 

Our results indicate that the effect of random component 
of input on random component of output does not last for a 
cast interval or more. This is indicated by Eqn. 4.72, for blasi 
furnace dynamics w|iich does not contain any term involving the 
operator B. This is of same order of magnitude as in the case 
of standardized input-output variables. However, the latter 
indicates the details of transformation and particularly the 
fact that the transformations involved delays on a part of the 
input, zero delay oh another part and transformations on the 
remainder due to time constants. Thus, for example, the outputs 
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silicon and sulphur contents of pig iron are influenced 
directly by the characteristics of the inputs in not only 
the cast interval just prior to that but also those of two 
earlier casts. The effect of inputs on hot metal temperature 
persists only in the current and next cast intervals. These 
results seem to be interesting. The order of significance of 
the inputs in the earlier cast intervals are also indicated in 
Table 4.5 and they are generally of the order of 12 to 35 per 
c ent. 

The feedback dynamics matrix, in both the cases, viz., 
in teimis of prewhitened input-output variables and in terms 
of actual input -output variables was found to be zero. This 
indicates that there is no feedback of the order of one cast 
interval or more from any of the output variables to the input 
variables. 


++++++++ 



CHAPTER 5 


SUIMARY, CONCLUSIONS AND SUGGESTIONS FOR FURTHER STUDY 

5.1 SUI-MARY A^TD CONCLUSIONS: 

In the present study mathematical models for blast 
furnace process have been reviewed. Mathematical models have 
been broadly classified as deterministic and probabilistic; 
the deterministic models were further subclassified as steady 
state and dynamic models, and the probabilistic models were 
subdivided into regression models, linear process models and 
time series models. 

A decoupled multivariate time series model has been 
developed for the blast furnace No.l of Bokaro Steel Limited, 
Bokaro Steel City, using data on three input variables, namely, 
sinter-to-coke ratio, blast humidit3r and blast flow rate; and 
three output variables, namely, temperature of hot metal and 
its silicon and sulphur content. 

All the variables , but blast flow rate, were found to be 
stationary. The blast flow rate series had two jumps and 
hence it was divided into 3 parts, each of which was stationary. 
All the variables were found to be free from cyclicity. In 
order to remove persistence (internal and/or external dependence 
a univariate model was fitted to each of the series. The 
resulting residuals which are normally distributed and serially 
independent are called as ’prewhitened series’. 
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A lEultivariate model was then fitted to the prewhitened 
series in terms of serially and mutually uncorrelated random 
components. The multivariate model is referenced to as a 
decoupled model since it decouples the within -the -series 
variability in terms of univariate models from the among-the- 
series variability in terms of the multivariate model. The 
advantages of decoupling are (i) as the parameters of univariate 
and multivariate models are estimated separately, their number 
is comparable to those of univariate and multivariate models 
and difficulties involved in simultaneously estimating all 
parameters are avoided; and (ii) the mo-st appropriate form and 
order of model can be selected independently for each variable 
and hence there is greater flexibility in the choice of the 
models. 

The method proposed for multivariate modelling is the 
generalization of the method developed earlier by Phadke et al. 
It uses G-ram-Schmidt procedure to obtain orthogonal vectors. 

A principal component model was derived in terms of prewhitened 
series and orthogonal vectors. A multivariate model was then 
developed by calculating the cross correlation functions 
between orthogonal vectors and multiple regression analysis. 

It was found that all orthogonal vectors form a family of 
serially and m.utually uncorrelated random components. Hence, 
the multivariate model was also a principal -component model. 

Prom the multivariate model, transfer function model 
between prewhitened input-output variables was developed. In 
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order to obtain the relationship between actual input-output 
variables the multivariate model was recoupled and then 
transfer function model was derived. It was found that the 
relationships betv/een inputs and outputs consist, of terms 
involving shift operator B which represents the time delay 
and other terms representing the time constants. The order 
of significance of the inputs in the earlier cast inteivals 
was found to be 12 to 35 per cent. The average time lag between 
input and output was found to be varying from 0.85 to 1.09 cast 
intervals. These results generally agree with the results 
proposed earlier by other workers. 

Peedback dynamics was found to be absent. This indicates 
that the control action has its affect atmost within one cast 
interval and there is no effect beyond that. 

5.2 SUGGESTIONS FOR FURTHER STUDY: 

In order to determine more accurately time constants and 
time delays, data are to be used at time intervals much less 
than one cast interval. Data on input variables were available 
at one hour interval. In this study they were averaged over 
the cast interval. For better results they should be sampled 
and used at shorter time intervals. Data on output variables, 
however, are generally available only at the time of cast. 

It is suggested that to notice the effect of operating 
variables on the pig iron process, one has to consider the 
extent of direct and indirect reduction processes which are 



functions of ore properties; its reactivity and the flow 
pattern of the gases and solids in the furnace. Therefore 
one of the important indicating variables about the state of 
the process majr be 'excess heat' above the normal requirement. 
This could be measured either in terms of top gas analysis 
(CO/CO2 ratio, E^/R^O ratio) or material and energy balance 
over different sections of the blast furnace. 

Furthermore there is a good correlation between the 
variation of the silicon content of molten metal with that 
of silica in the burden. Using these and other pertinent 
variables, more comprehensive models for system dynamics of the 
blast furnace process may be developed using the decoupled 
multivariate modelling approach. They may then be used in 
developing control models and control systems. 
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APPENDIX A 


MARQUARDT ALGORITM FOR NON-IINEAR LEAST SQUARE ESTIMATION 


Let 


^ =( € , ^c;)denote the parameters of the model, 

that is = ((p^,q) 2 » • • • »9p» ‘ ‘ * ’^q.^ where K = p+l. To 

start with the initial approximation a convergence para- 


meter e and the parameters n and P 2 are specified which const^’ain 
the search. During the search, the values £(t) and the derivative 




3s (t) 


?£ 

a 'ji 


are evaluated, niimerically, at each stage of iteration. The 
derivatives are obtained from 


X 


^ ^ '=^ 

»••••» 'ott- n -' 


i,t “ "^IjO’**’’ ^i,0 




1-1 p. , . • . , p, ...., 




' 1,0 


i,0 


■K,0 


)]/ 6 . 


Stage 1 ; With e(t) and x. , supplied from the current para- 

1,1: 


meter values the following quantities are formed: 


1. The matrix 

i = [ hi 3 


N 


where 


A. . = ' 

1 ' 


. = X . . X . , 

J > i,t 3.t 


t=Q 


2. The vector g with elements g^,g 2 » 

N 


% = 


X . e(t) 


t=Q 


where 
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3. The scaling quantities 


Stage 2 ; The modified (scaled and constrained) linearized 
equations 





are constructed according to 


A. . = A. ./D.D . 
rO 13 1 3 


i 3 


A. . 

11 


1 + Ti: 


gf = gi/\ 

The equations are solved for h which is scaled back to give 
parameter correction h., where 

h. = h./D. 

3 3 3 

The parameter values are constmcted from 

I = + h 

and the sum of squares of residuals is calculated using, 


S(* ) = [£(t)] [e,(t)f 

Stage 3 s (i) If S( -^ ) <3^ S( the parameter corrections h 

are tested. If all are smaller than e, convergence is assumed 
and the matrix A~^ is used to calculate the variance-covariance 
matrix; otherwise is reset to the value x is reduced 

by a factor P 2 and computation returns to Stage 1. 
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(ii) If S(?- ) S( s ), the constraint parameter ti is 
-f '*0 

increased by a factor computation res-umed at Stage 2. 

An upper bound is placed on n, and if this bound is exceeded, 
the search is terminated. 

Ivhen convergence has occurred the residual variance and 
the variance-covariance matrix of estimates are calculated. 
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APPENDIX B 

GRAM-SCHMIDT ORTHOGONALIZATION PROCESS 


The process generates a set of K orthogonal vectors 
j 3 ^(t), k=l,....,K from a set of K linearly independent vectors 
£ (t), k=l5....,K by forming linear combination of e (t). let 
a^(t), k=l,....,IC be a set of orthonormal vectors, then the 
orthonormal set a(t) has a property that 


(a^ (t ) 'cc^Ct) ) = 1 j = k 

= 0 j k 

the paranthesis (•) means the inner product of vectors a.(t) 

3 

and a^(t ) , 

(a^ (t) va^^Ct )) = ^ “k^^^ 

i=i 

To start x^ith Gram-Schmidt ortho go nalizat ion process, set 


P^(t) = £^{t) 


j|Pl(t)| I 
a^(t) 


= (P 2 _(t) -p^Ct) ) 

_ Pl(t) 


X 

2 


The second vector unnormalized orthogonal 

vectors is found by a linear combination of a^(t) and 
that is, 


P2(t) = ejCt) - 63_2 a^(t) 


where 6^2 is a constant. 
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Forming the inner product with a^(t) we find that 

(P2(t) )) = (£2(‘t) ^ 

Since ^2^'^^ p2_(‘t) 3 .re required to be orthogonal it is 

necessary that 

(p2('t) ' ^ 

By orthonormal property 

(a^(t ) • a^(t) ) = 1.0 

Hence 

6^2 = -oc^Ct) ) 

The third vector p^(t) in the set of unnormalized orthogonal 
vectors is formed as a linear combination of a^(t), a2(t) and 
a^(t) expressed as 

p^(t) = e^(t) - 6^^ a^(t) - 62^ 

By taking inner products 

(p^(t ) -a^Ct) ) = (e^(t) -a^Ct) ) - 6^^(a^(t) -a^Ct ) ) 

-62^(a2('t) ’cc-LCt) ) 

(p^(t ) •a2(t) ) = (e^(t) •a2(t)) - 6^^(a^(t) •a2(t) ) 

Since p^(t) is required to be orthogonal to p^(t) and ^2^"^^ 
have 

(p^(t) *a^(t )) = 0 and (p^( t) •a2 (t ) ) = 0 


By ortho normal property 
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(a2(t ) •a2(t) ) = 1 (a^( t ) -a^(t ) ) = 1 

and (a^(t ) •a 2 (t ) ) = 0 
Hence 

*^13 ~ (o>:3_('t) ’s^Ct) ) 

*^23 " ‘EjCt) ) 

The general expression for calculating p.(t) from £.(t) becomes 

3 3 

3-1^ 

pj(t) = £j(t) - (a^(t)*s^(t)) a^(t) 

i^i' 

or p.{t) = o.(t) - ^ 6 ^^ a.(t) 

i^ 


where 


a^(t) 








3jt« _:ic*,jjc.,. jj<*= ^»=s|e 


'J.P''( hC!X 


s^i.v;, I*-- ,-s. ;.J, j;:. sv. :}£•*- *-»!5 s-=- :|e. sjc^jjc., lij!... s}!,,, aft* >jt^ «, :J5„ „ 4: i , 
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■' P P E N D ! X D 

»-. t . r* « « «• t «a(t.3Si.r€ «•..’«.■?? «• frs«-’« * iKS«,i-ir»fc.*»»:- 

f f f-’ P U”^ ” P FRfGHiaMS *’5^^=4<=»* . : 


:^. t >:;!’• y- ;*■ ^ 1 ' 5:-^ 3t|« ^ # :f« 5^ ^ sfc # ^ ^ ^ ^ # sfs ^ ^ ^ ^ 

••f'.:.f-'f' i r> i ir CaLCULATF AUTOCORR flat ion and partial AUTOCORRPL/Ti dk 

FUNCTION OF THE SERIES «* 

I i; 5?c ^ ^ ^ ^ ^ ^ 

?:-■ J ^ *"T nf.ic , 5?5ci?^e :^ . 

'fvFfT) AUTOCOVARIANCE FUNCTUON*. 

'fFin oc AUTfX.CRRFLATinN FUNCTION® 

o/LFr,n *». PARTIAL AUTOCORRELATION FUNCTION. 

'■' KUNFFR OF PROBLEMS. 

’ roo r. UPB'-F OF DATA POINTS. : • 

ore ^FAN OF THC SERIFS. ■ 

"V ,.c VAPIANCFCOF THE SFRJES. 

n ... STANDARD D^VIATIDN OF THE SERIES, 

'U>/4MI) ... STAFXARDIZHD SFRIES. V ' 

••'! ’r .„c CK^RFF OF DIFFf RENCING. 

^1' {eJ:, vx 55? * ,3^ ♦ # ^ ^ ^ ^ # # jIc Ic ^ # 3^ # ^ # # 

C'JFtfPJnN WWtlOS^i-) 

n URL^ PRCCISION WBAR(1^5^ ) 'P 

ri LRIX PR! Cl SION Wdrsr ),ACVF< 1’3M,ACF{1DC) tPACFiSOfSO) 
nrUPU* PP' Cl .'JON SBARf SMBAR,SSBAR,SVBAR,SDBAR 
r( ijBL'- PP’ Cl SION SM,SMF AN,SS,SV5SD,SNUM,SDEN 

' CP'='^'' 

! prr, F=f'f'> ' ■ ' 

n AD 3; I»MT ' ,'' ■ ,: : ■ ' 

'■'C II’ • ri7S=ltMT ' ’ .' 

ri / f)3i ! »t-' 

AO Jf 2, {WWU ) ,1 = 3 ,N) ■": •- ' 

DC I I» 1 ,N’"’ ” ' ■ ^ 

1 W{|)«D 8 LF (WW!I) ) 

pp I f-;Tl (:?> 

IF{NTS,Fo.1) print 121 ' ^ Xp 

"F{N^F.PQ.2) PRINT 12 ?. '. ^ ; r. 

IFCNTScfO.?) PRJNT 123 ’ ^ 4 ;;::;/ 

^FINTF.e Oo-f,) PRINT 124 

I Ftf!’'? .FQ.T ) PRINT 125 • ,-, 

J f ( M* 5 p ! 0 , e- ) PR I NT IZe , , . '"U ; 

CTLCUl/TI’ WFAN AND STANDARD DEVIATION OF THE ORIGINAL SERIES. " ^ IrP/;’ r’' 

'■M = ' .'“T''’ : 

' "r.l, f, 

'‘ ^f'i = ?F+wn ) ■ ' . ;■ :' 

'■f,’’.' N= 5 t'i/UBLF IFiCAT(K) 3 , 

“S = 1 ' « M.” ' ■ ■ .; . . ■ . 



^ " 5 = c5 + ( Wf r )• SfT ‘ ■ 

-- f FLr.AT (N) ) 

' i - • "C! 7'f-V) 

r p r-^. . , 

■■ i ] I” ” ] (* i . , , _ : ' 

zri ,MlfS.VfSD 

C ■ !^ ’ SfRJIF E^Y SUBTRACTING MEAN AND DIVIDING BY STANDARD Di' VI i. TI CRi, 

; Mf / I (T ) =-c^ (T 

■ r,f : 5? , . 

!• 3 

^Ffr.TiFcroci )Gn ir 

.1! I".;:,,. ’!■■■•'■■,, ■ ^ ■ ■ 

? '^F / r f T ) =wr.AP. (I + l )- WBAR{ I) . 

/ 7! r- ' : f iir 

C CntUL/'"r Ml AN AND STANDARD DEVIATION OF TRANSF3RMFD AND DIFFEREMCFO, IF 
C (if-JF ofiio? ADC ?), STRIESa ^ ^ 

^F/.P-f D( 

'-.f, t'f J=1 jN 
;■ ■ '•rAF=$P.AF:+WBAP( 1) 

' M- / R=SEAF/DBLE{ FLOaKM ) 

"SFAF^i *t DP 

'1 I=lfN ■ 

?.l S5PAr = SSBAR+ (WBARCI J-SMBARIYYZ 
'' VfUP=SSBAR/DBLS IFLGAT IN) ) 

^r)P/P.=DSQFT{SVBAR) ' 

c caiculatf autc correlation function^ , ' f 

Dr 6 1*3., 3 AC F ’ 

JcvF{i)=' *on^ ■■ ' 

f P, ={.,!. J . : : 

nr *? J=3 .NM ' . ' . ■ 

J J = J+I ■ / ■ ' 

" ' CVF(J)=^ACVF(I) + (WBAR{ J)'»SMBAR)>«'{WBAR{JJI‘“SHBAR) 

«CVF {n=ACVF (n/ORLLCFLOATIN)} 

• CF { ! ) “AC VF( I )/SVBAR 

- ' ’'.V':-'’''. '■ ' 

i'f 5}; * :}•. Jt: ii: ; 5 s : 4 e * :^ * :*; zje 3j: 34: # ♦ :4 p # Tjc ;{c ♦ * * # ♦ s}t ; '#( 

C r M- DLL P/RTIAL AUTO CORRELATION FUNCTION , "/r 

^ 3,v s!< « 5?£ ♦ * # 3js * # 3J( * 5j! 5{( * 3jt # ♦ # # !je 3jt 3#: ♦ sic * ifc , ■ j_ lj|«^ ' 

p;rF(i^,i)=ACF(i) ■ , ^ 

PC n l=2t ipacf 

l L-‘f ' 



' J-L- J , : - 

Mj.= L. 1 

■ t i-^i-ur', + P. CF (LLLr J)*ACF{LJ 5 
o ri + P,M. FCLLL, J)=<'ACF{ J) 

^ Fd; L) = {^(:F(L)-SNUM)/( 1 ,DD 0 «»SD 6 N 3 
" ; It J=^ 1 , ,LL 
IJ J = L‘ J 

1 -^TF (Lf J)=PACFCLr 4 sJ)”PACF(L,L)’S'PACF(LNtLJJ) 

] 1 Fi r.7 -r )!'■ . ' . 

Of, -■ A -1 ^f:nj F 

Of-- |i;i" ; S'P : ' . : 

•"nr," 1 F 9 . fPBAR, SVBARjSDBAR : 

C PP^n-.-'-H'' dTC COVAPIAFJCE FUNCTION) 

PFim itt'- ,t , 't', 

■ 21 .) = ! ,lACFt 5 i.' '■ '■ t ; . -VF 

’ ; = j+.',q • ■ , 

?t H'n (ACVFn)» I = J, 1 F) ■ 

C PRM'T TU’ AUTC CCRFFLATION FUNCTION® ■ ■ ' 

i’p.i fr 1 52 ■ : r., 

"•n 21 J = ). ,IACF,&C ■ 

Tr=j +*.9 . 

2 ’ PPdd ) 5'1 ? CACFd ),I=J,IE) 

C PFIIT "HF PARTIAL AUTO CORRELATION FUNCTION® 

, ' , ■ . ■ ■ ■ . A;'. _ A • AAV" lF-A:v;.;'v-::,viA!^ 

PRINT 1 - 5 A ■ ; ^ A ' n: A^ A::nAA* 

no 1 ^ J= 1 ,IPACF ,50 . ■ : A,. . :P: ,AnAAnsA||^^ 

n =J+A 9 ■ 

4 A rPTf'JT m ? (PACFI I, Hf I=J»IF) ^ A 

32 rONTTNUF , '■ ' ' '• > A A' 

3 - AH' Cf.'KTiK'U'- A ,A- , A;,', ' • 'A;-; A'i 

Jftpir » OUTPUT FORNAT STATEMENTS® '' ’ , 

Sli :(• 5'r A * v « # 3^ ^ * Sjt ^5= * * :«£ !$: :{c * ^ * :je :4i; ^ jjt ijt + ^ !*: j}c jjf ' -t J 

/rp tjf FA Cirm' IJMT JOB irRMINATFD ■ •, , , y; ' .jy 


i*i" A. 
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-'fu. / TO CALCULATE POWER SPECTRA OF THE SFRIFS »» 
r. SP C"F';’L A,f!ALYSlS OF TIME SERIES ... 

r -'-HV' prnGp-, 5 f- uSf'S hA MM! NG-TUKF Y METHOD OF SMOOTHING. 

C th' T PFOr-HAH ALSO CALCULATES CROSS-CORRELATION AND AUT0-C0KRELATI3N FUMCTIO? 

sV: :y. f s!- s!< St' # !^r * 4: * # ^ # sjt 3!' 5(5 ♦ ^ ^ * * * # * * * * 

s-j r-'0*'A~7rT.IS, 

wp(n autcovarj AHCO function. 

?F'{J) AUTCCORRFL&TION FUNCTION. 

FLP{ ! ) ■ RAW SPECTRA, 

UPC I) smooth SPECTRA. 

- syr. li-.er, sy -= -^<iT syo-.- A. .*«!»#.. s!<», 5§:«»>5{ec-. 45e~'3t5«^ #••• 5};«» I. :(: >.:((» 5je..:{i:««. 

n I Mr NSI or.! X ( 6 t35f)) ,TPO50 ) , SP { 353 ) tFP{ 350 ) , GP t 350 ) , Z ( 6, 350 ) » 

’ 7Z(<' ) ,CP (A'O) tWPCA!^) ,RP{AU),UP(43 ),FLP{40) 

■ ;-n rifM^'jNVAR 
Pr 009 LM=1,NVAR 

. TF {LMoLF.?) read 102f( ZCLM,Uf I = l,NN) 

" F C LM.'*Q«3) GO TO 1111 

TF CLfUG'".^!-) READ l'2i( Z ( LM, U , I = 1, NN) 

90 TC 9^9 

1311 FfAD 1 2,{ Z CLMj I) , i = l,4D ) ; 

T"AD 104,( ZCLM,I),I=41,A3} 

5'/.D 1' 2»{ Z <LM,I ) ,I=44,NM) 

f,> 9 rcHTir-tur 

>FAD 1* IjNPROB ■■ ■ . 

nr 600' LNM=1,NPR0B 

np :^'"0 17=1, NVAR ^ ' 

OP IY = 1 ,NVAR - . , ,,J- : ‘ 

SUNX=( .P " ■ O'O- 

suMY^o.u ' ' '-'O' 

3|r Iff* »* f** ^ ^ 5^5 *•' J§Cw3 5 ^s*' 9 *§C*»3^«f3 

.C C^LCULATf MOAK AND STANDARD DEVIATION OF THE ' SERIES. 0, ^ 

DO 3 r? i = i ,NSAMP 0 ..-P'.' . ' 

SUMX = SUMy. + Z<lX,I ) ■ ‘ ' 

FUMY=SUf-!Y+Z(lY,n 

' r.f 'Mjrjuf' 

’•MrAP=SyM:*./FLnATCNSAMP) 

YPf AN=SUHY/FLnAT{NSAMP) 

" v''=i o! :■■■ * ,: 

C VY = * rj- - ■ - 

PC 'U 7 J=I,NSAMP 
" VX=SVX+C ZCIX,! l‘-XMLAN)=^'=i'2 ' 

rvY = SVY+< Zf IYtIJ* YME- AN)^*2 
: r -(^NTTPui . : 

5 v3'=sv! / floatcnsamp) -’■pr 

f VY=r VY/FLOATC NSAMPI 

5r yx^oRTC svx ) 


Y=rGr‘" ( 5vy ) 

:- «• fYYI Y^'FAN,SVK,SVY,S0X*SDY 

r f < pG ■■. 7j tk: sopi fs. • 

5<;.. 5^ - J' •. i’ r- 5*‘t ,*’t,f s?.f x^,^r -^fw-::kb-:^r’'f:t 's^m:^m':$:mm:^m :^m:^m» :^m :^m ^ :^m m m b' ^rr 4^ »* Jfc.* 

nr '." T=], ,Kfj.|yip 

■ f’. -n=^fZ(l>»I)'YMFAN)/SDX 
' t‘V,n-(Z{TY,J)*-YM»AN)/SOY 

; ' ■ :i ! ’'’ff'j! ' 

"1../ = PL/' G- 1 

sj? S,? :Jr A 3?*: ^f yc 4 5^ * ^ 3^ ^ ^ 5^c ^ ^ 

CZLCnL;-^r TH?" ALTOCOVARIANCF AND AUTOCORRELATION FUNCTIONS* 

3% A s'" >!"> ri*: jte # # 3^ # ❖ . 3|c ^ ^ ^ ^ ^ ^ ^ ^ ^ ♦ 

■"PCI >=r 
npH )='‘,^ 

TFC] ) 

^P<:< 5=' ei 

DO Jf. 1=1 sNSAMP 

rp(] )=SP(i)+y{ix,i )=4'*2 
''■p{l)==Tp{i) + y(iK,T ) 

":pci )=GP{ l) + :v{i Y,i )*=fr2 . ' 

FP{i)=FP(l)+YCIY,T ) 

1 .'■« PONTJ^jll!: ' ■ 

mi=mlAG+1 

Of 2f i:=? ,M1 

■ . 

K = NSAMP-’I+2 

Tp(i ),TP( J)' X(IX,J) 

^P(I )=SPn"*X )-X(IYt J)=«‘=«'2 
FP(I)=FP( J)-X{I Y,K) 
r,P{I )=GP(I-1 )'»X{IY,K)^’S'2 
?•■ CONTINUi: ' 

MLAr,l=MLAC-+l 
Pf' PI" 1=1 tMLAGI 
MMINP=HSAHP- I + l 
CP(I)=taU 

on 26 J=1?NFINP : ;.■ 

Kt =J + I*’l : 

PA rp {n,=CP(I) + (X(IX,K2)=«'X{IY, J)5 
WP(l)=CP(n /FLOAKNMINP) 

P[',-UM=FLOi TCNHir4P }=<<€? (I )“FP( I )*TP( I ) 

PlVf.'l =SQR‘^ C (FLnAT{NH!NP)*GP(I) J^FPd )=*'#2} 

- iGhipsFQp-^C (FL0AT(NMJMP)!^SP( I) )-*TP(I )**2} 

^FCJ f ^‘’,IY)=PMUM/«RDFN1=^RDFN2) 

C CAICSJLAAI ""HG PAW SPFCTRA. . ■’ , 

Pi, '’ 9 ) =1 sML/GI ' ■ .V' ■ : ■ ■ 

PI PC T )"' . . ■: : ■;. , 



T'' j=r 

• K!.f>!:')=FlPCl)+2«=4'WP(J)#C0S{3,l415927>«'FL0ATUI-l)*{ J-l))/FLOATCMLAG 

■J ) ) 

-■ (’> = rLMI) + WP(l}+WP(MLAGl)’«'C0S(3»14l5927=!=FL0AT(I-l)} 

C 5! -’H!' 5 ?f CTRJ?, - 


M;- fi)=' *FLP{2 )+r’c?4=!‘FLP{l ) 

{wL;-f-/i )=nc^-6’^FLP{MLAG)+0t,fj4=«'FLP(MLA61F 
"^r I=2TML/iG . 

MF.n)= *?G*FLP( I- 1 )+C.54*FLP (I )+&*23*FLP( I + l) 

nc -’f . . 

T p = }. ^ 

7'- Ppjri-"' l(",IP,WP(I5,FLP{I),UP(nfRP(n 

:•■.( I rr f'T-’-f'uF 

nrr rf.F’-rTi-;ui' 

IM F( fPF19.8,F16.8) 

in FCiFMAT U'' If ) 

1-'? FnPP.AT(!jri6»7) 

If:.. rFFM/'^{'’n6.^) 


?' / n PMf T(1!",*XRrAN=’«',F2P,8»5X,=!'YMLAM=«‘,t2!V*8,5Xi=«'SVK=*«'»F2f#,8/flX,»SV 
’'Y=*,! 2('.r ,'"X,^SDX=*,F2C,8,FX,=«'SDY=*, F2D.8) 

9f-. Ff'FMAT(2Xf* SPFCTBAL ANALYSIS OF*, 14 , *VS*f 14, //’i' P 
]CrVAP. PAW SPFCTRA SPQOTH SPECTRA *//) 

5Trp ■ ■ . ^ ^ 

'TT,, 



5^.. MI'S, 


■ , 2 'j ■' C 


s^esfejji 


OBJfcCf PRGGR^F IS BEING ENTKRfC INTO' STORAGE AT 05 HPS, 




PROGRAM NO. 3 TO CALCULATE THE INITIAL ESTIMATES DF THf FARft^■i"Tt• PS cn 

3^m* 3?c«ro 5^«a ::^«» 3$c«®i 5fc«» 5^« 3^«» 'sfScw 5^* wt *» m 

*=«f* this program calculates initial estimates for univariate stochastic 

MCDFL of a time SFRIES. THE ALGORITHM IS GIVEN IN P.rFFRlNCt NO. o5- 

m, :^tm S^easSjS to t^m'z 5^« :^es 3^rs #c ■? iSs*» 3jC«» 3|etes 3§S»te 3§;«, 3^«» :|i:*« ^ ^ sea :§C t. # « # t ^wr w* 

NOTATIONSo.o.o. 

PHKI) ... INITIAL ESTIMATES OF AR PARAMETERS. 

*** THFTAd) ... initial ESTIMATES OF MA PARAMETERS. 

ITMAIC ... MAXIMUM NO. OF ITERATIONS. 

*:«cs*c NPROB . 0 ® NC« OF PRCBLEMS. 

#=«'* SM ... MEAN OF THE DATA SERIES, 

*** Cd) ... AUTCCCVARIANCES OF THE SERIES. 

*=«'* MP ... ORDER OF AUTCREGRFSSS ION® 

*** MQ .®® ORDER OF MCVING AVERAGE. 

*=«=* NDIF ... DFGREF OF DIFFERENCING. 

3jc&)5c ^^»5^c«sa4c«»3§:«^ :^»us^r* . ::^«a ^«»3§cr-i 5§f«:^ws #«■ ^ * S 

C UNIVARIATE STOCHASTIC MODEL PRELIMINARY ESTIMATION (USPEI 

DIMENSION PHI (6) ,THFTA (6 ) ,C( 10 ) ,TAUI 6) t F( 6 ) ,X { 6 ) *H < 6 ) ,CH « Ifi ) 

DIMENSION A (6,6) ,T( 6,6) ,T1 (6,6 ),T2{6f 6) 

FTA=C*061 , ; ■ 

ITMAX=1(1D ■ , - ' ■: 

READ l,NPROe 

DO 600G L=1,NPR0B ' ■ 

READ 2,NP,NDIF,NG 
READ 3-,SM 
NPQ=NP+NQ+1 
READ 4, (C (I) ,I = 1,NPQ) 

c print the input 

PRINT 204 ,NP,NDIF,NG 

PRINT 205, SM 

PRINT 2C6, (C (I),I=1,NPC) 

IF(NG.GT.C ) PRINT 4C:l 
IF{NG®GT®0) print 149 
NQC=NQ+l 
NPP=NP+1 

IF(NP.EQ.€) GO. TC 3L0 

C CALCULATE THE INITIAL ESTIMATES PHI OF AUTOREGRESSIVE PARAMETERS. 

DO 6 I=1,NP 
IQ=NG + I+1. 

X(I)=C(IQ) 

DO 6 J=1,NP ■ 

IJG=IABS{NG+I- J)+l 

6 A( I, J)=C(I JG ) : ^ 

G SOLVF THF SET CF NP LINEAR EQUATIONS. . 

C SUBRCUTINF. MATINV calculates the ' roots OF THE ,SIMULTANIpUS ‘ equations# 



CALL FATINV {A, M F, X, 1 , DFTER } 

DO 15 I=1,NP 
15 PHI{I+1)=X{I> 

DC 33 J=1,MCC 
CH<j)=ron 
DO ZA .I=1,NPP 
DO 24 K=1,NPP 
I JK=I+J“K 

TFUJKoLForj) IJK=IABS( IJK)+2 
24 CH{ JJ=CH{ J) + PHI (D^PHICKl^CCUK) 

33 CONTINUE 
GO TC lOf? 
im DO 42 J=1 ,NGQ 

C CALCULATO INITIAL ESTIMATES THETA OF MOVING AVERAGE PARAMETERS® 
C THIS UTILIZES NEWTCN-RAPHSON ALGORITHM® 

42 CH(J)=C(J) 

leO IFING.EQ.O) GO TC 2400 
C THE ITERATIONS BEGIN HERE® , ' 

IT=1 

IF{CH{l),LT.e. ) GO TO 111 

> GO TO 112 ' 

111 PRINT 222,CHC1) 

GO TO 6C0r- 

112 TAU{1)=SQRT(CH(1)) 

DO 51 1=1 ,NG 

51 TAU( I+l )=C®0 
240 DO 60 J=l,NCQ 
F( J) =0,0 
JQ=NGQ+1~J 
DO 69 I=lf JQ 
I J = I +J’“1 

69 F(J)=F( J)+TAU(I)*TAUIIJ) 

F( J)=F( J)* CH(J) 

60 CONTINUE 

C FORM THE MATRIX T CF TAUd). 

DO 78 I=1,NCQ 
DO 78 J=I,NGC 

Ti(I,J)=t}.0 ' , - • 

78 T2(I , J)=o®0 , : ■ ■■ ; 

NNQ=NCQ . ^ 

DO 87 I=1,NQG ^ 

DO 96 J=l »NNQ . 

11=11+1 

96 T1{I,J)=TAU(II) 



NNQ=NNQ-1 

II=I 

87 CONTINUE: 

NNQ=NCQ 

1 1 =c:; ... 

DC ICf- J = I,NGO 
DO 114 1=1, KNQ 
11=11+1 

114 T2(J,II )=TAy a ) 

NNQ=NNQ-1 
I I=J 

IDS CONTINUE ■ 

DO 123 1=1, NCQ 
DO 123 J=1,NGQ 
123 T{I,J)=Tl(I,J)+T2tI,J) 

C SOLVE THE SET CF (NQ+l) LINEAR EQUATIONS, 

CALL MATINV (T,NGQ, F,l, DETER ) 

DO 132 1 = 1, NCQ 
132 H(I)=F(I) 

C CALCULATF MOVING AVERAGE PARAMETERSA 
DO 141 J=1,NQ 

141 THETA(J+1)=‘ TAU( J+D/TAUCl) 

PRINT 150, IT»(THETA{J},J=2,NQQJ 
DO 159 1=1, NCQ 
159 TAU( I)=TAU{I 
IT=IT+1 

IFIIToLTelTMAX) GC TO 160 
PRINT 600 

60 TC 240f; ■ ^ 

C CONVERGENCE TESTING® 

16C DO 169 1 = 1, NCQ ' _ , ' / . {« 

169 I FIABSIFCI ) ) .GT. ETA) GC To 240 . , . 

C IF ALL Fdl.LT. ETA ASSUME CONVERGENCE. ' ‘ . 

C CALCULATE; INITIAL ESTIMATE OF OVERALL CONSTANT rHETAIf ). ; 

24C:fi I FtNP.CG.C ) GO TO 330 

PHIS=C.D ■ ■* 

DO 177 1 = 1, NP 
177 PHIS=PHIS + PHHI + 1) 

THeTC=SM=«' 1 1.0-PHIS) 

GO TC 2Ur 

330 THFTC=SM ^ 

■■ ■ ■ 

C CALCULATI' initial ESTIMATE OF WHITE NOISE VARIA^fff 



I F(NQ.GT./',) GO TC 186 : 

PHIC=f.6 
DC 195 

195 PHIC = PHIC + PHI(I + 1):S=C{I + 1) 

VAR=C (1 )-PHIC 
GC TC Bie 

186 VA.R=TAU {!)*=» 2 ■ " ' • : 

31f- IFCNP.CQ.G) GO TC 21P 

PRINT 2U8,{FHni+l)tI = ltNP) 

21P COf^TINUF 

PRINT 211,THFT0 
■PRINT 212, VAR 
6e00 CONTINUE 

3^ i» jjco ,5 |:«w :§c*» sSs«» 3|c*» 3|s» sJc w» :fcar- 3|c«. « :§e ^ • :^w ^ 

FORRAT FOR INPUT - OUTPUT. 

1 F0RMAT(lf)I5 ) 

2 F0RHATC3I2) 

3 F0RMAT(5F12o7) 

4 FORMATC 6F12.?) 

149 FORMAT ( //5X,*ITE; RATI ON*, 12X,*THETA{ I )*// ) 

150 F0RMAT{5X,I5,5E2067) 

204 F0RMAT(//1X, BORDER CF THE PROCESS IS ( *,312,*)*) 

205 FORMAT! //IX, *THE MEAN OF THE SERIES IS=*,F12.7) 

206 FORMAT! //IX, *AUTC CCVARIANCES OF THE SERIES ARE*,5F15«7) 

35ACH ITERATION ARE AS FOLLOWS*//) 

208 F0RMAT{//1X,*INITIAL ESTIMATES OF AUTOREGRESSIVE PARAMETERS ARE*, 
1E20.7) . , , 

211 FORMAT! //IX, *INITIAL ESTIMATE OF OVERALL CONSTANT TERM IS * ,EIS# ?| : f’r • 

212 F0RMAT(//1X,*INITIAL ESTIMATE OF WHITE NOISE VARIANCE IS *,F15.8) 

222 FORMAT! IK, *TAUO IS NEGATIVE TAKE NEXT PROBLEM*, E17 .8 ) < ' 

401 F0RMAT!//1X,*INIT|AL ESTIMATES OF MOVING. AVERAGE PARAMETERS AT 

600 FORMAT{//5X,*NOo CF ITERATIONS EXCEEDS THE LIMIT HfNCI PROGRAM, I S , 

eTERMlNATED* ) ■ . ■ 




STOR 

FND 





" ''f, 4 FINAL FSTIMA.TICN OF THC P/iRAMFTL'RS * 

:■ FTP. ^-^'FAHrTcR rSTlRAflDN 

I. • : : 7 - algorithm for non«l inear least squarfs 


d/*> #is3|;e 3^ » # »» 3^ » 5^ ^ ^ ^ ^ it i ^ X # ■■ 4s- 

C If FRf GRAF # ♦ 

s/f r ’t%- :d:.m id* ::^a» 3§£*r' ::^»9 3|f,w3^ i» 3^*' :#‘i 

>,-. i '• T / - I C f-' S 

f’F f''!-!! HR RCVilC- AVIP/:G0 PROCtSS. • 

F.r r j. r’ p (4 / URRRf CP« f ^lOK « 

L l: m‘ f.p th:., sf rirso 

irjr f'f-r'i [T u. FFRrr I Cl NG, 

> f 3 :;' J! ■ <'.! F VfiRJ f Bl' . 

KK jf ■ ' iael: 

KK=. ■’ PFOGR6F iS USED BY MAIN PROGRAM TO CALCULATE AACT) AND AAJl) 
SGUAR* PHI AND 1 HcT A E^JTERrQ AS ARRAY BETA. 

KK = Z PFOGRAM is used BY MAIN PROGRAM TO CALCULATE AAdl FOR 
p, TA crifTURBrD BY STALL QUANTITY DEL* 

'iH:'', H) /ROLY OF FtVIhG AVERAGE PARAMETERS. 

PHiU) /,RF'/Y OF AUTOR GRESSIVE PARAMETERS. 

AM") ill-': VARJ./ TL RSSIOUAL SFPIFS. 

05.' r.tlBL, PR:;CISiCN variable , UNCCBOINITIONAL SUM OF SQUARES® 

j; « ^ f » 7*C, t » ^6f'* Jlfwt s^irt 3^#» 5 I|£wp «i 5^ A S^i» #''W» # j 

PPi CJ SION USSfUSSO ' ^ 

■ ( IIP!/ PPM I SiO(\ URX 

• C’ F r.pTGiN ■ - ' ' 

' f"" F M /'CF (6' ) ■ 

" v.|;M'L WVMSTMfP (575) , WOTS) ,Z(375),PHl{4>,THfcTA(4>rBETAf4) 

'M'M; RSIUH WWW(375) . ... 

r M'MfJDN WCi:75) 

rM'M.c.HjH ::>, M,?TS),A,AO75),GC4),BB(315I,SC4),BETA0C4),B(4ti) 

. r: MS I'SMiC AAA(4,OtAAAS (4,4) ,DI4),GSC4>fR(4#4),V{4t4)fH{ 4>,HS{4I 

M PMCr' wu 

r /. n : ,nppc8 

:t : f' Lf^ = l,f'PRCF 

' I I- iJ ] 1 , I’-J : . .. . . ■ , •:' 

'■".=R ■ ' ■ 

■ r. T ?, ( Z(I ) ,I=1,L) 

C . MO!./-, p; At. AMO STANDARD DEVIATION 

'■■ ! R =: F „ 

!■ . ^ '■ I = : ,0 ■ . ': 

I ' '■ M r S'.i;'i4 z ( I ) 

' i M- .r-.Mir'/FLOAK A) ■ '■ " /'■- 

.' .:.■ ■ ' ■ ■' ■ .. ..' : 


l=! ,N 



* ' ■ . . ^ — r (Zl I )• 

V’ /f uf'f T(M1 

f; cf- "( ' \//.) 

' ’ ' " = l,f 

; ' a f T )' CM/SD 

5 ' ■ j-’rF.u 

:■ /'■ ’ yi.,P f 4, 

■^r'D,r,To' ) K AD ■ijiPHKI ),I = 1,NP) 

• '! rr-]-. ,'■) 

•FD’C.^T* ) F AC {^H: TAJT) ,I=1,NQ) 

' F‘F4 "? '^Fl' «7) 

"F^'n’T'D F 

i] . "1 F : 7 /-I f FK 

i'. : f r ) - wwiv; ( n 

Tr 

F .' 1 . ■■ 

i j ' 4' r ( :! ! = w u W I I ) 

■> . y| 

I=I 

1, - {T (T+i)» WCd) 

. ?. ( f . I -! jf'N , 

11/ Wf { ■ )=WW4'(1 ) 

''■-■■•f'T • ? 

7. 'FF !=I,^ . ■ 

■ >H (1 )=Wini42 )-2,y*wr:(J + l)+WD(n 
i: MU' ■ 

. V~ o 

of 

. ;:-t ] i=i ,K 

"i- ff =SF+W(:>{i ) 

"•f sSF/FLfAT (K) 

i ( ■7'Tl T«: ,M 
•71; 7v= 5 v+f wrii 

' V^f- V/Fl.fy-'" (fU 
"r; = 7CFT!SV! 

?!!7 ■ sSfSFV 

T pf'. [ :F .r C. ) GO "^0 ?D<'8 . 

nr ■■ r / 1 = 3 ,1,1 

■ ‘ ■ ’.i t ; ) = »W1 ?l )• SP) /sr 

prf:'' I- *t 
.'■ V!r' 1"'=' » 
rc ■ I =1 ?N 
:■. " 1 r; 7, = r. FP/ k + wn ) 

n i- t :-=-r: FD/' ?./ FLOAT (N) , 

nj ' , i j = 1 , |\, ■ ■ 

1< : VP r !' -T'. VP.AF.+ {W ( !)•■ SF'BAP 3’S'*2 

'■ VP' VL.‘.f/FLCFTfK) 
n? : > •" '* ,PFBAR ,PVtAF 



- 1 t . . r - ! Si:RIE:S==«<,E2r oT/1X,=«'V/RI/.?^jC; 

. ) 

■■ t ■ r ' ■» 

' : • • j =1 f r> 

= 1 ,; c 1 | . ■ ■ ' 

■ : ' ’I ''ij^' 

' ‘ ‘ IE 

' ■ V;v, ! T ) r , = ' ..■ 

^3 , :■ ,, ■ ■ : 


" i I ‘ ' ■ , . : 

I inr^ 

- 

C tjP n- TA 

’ r- {-Pc' Cot ! GCi TC 6 1 
" I ■' T - itl-.P : 

i . [' -'f- {1 3 = ^HI {7 ) ... 

6 i " Q,? ) Gli TC t’ 3 ' ,: , : 

T ,f.C' , . . ■ 

. P'' = r P+i ■ ■ , ■ . 

f,; , r,! T/, (cp, }=THcTA{J) ' ■ ■' ■ 

/ ■ ’ .i-('P + .iC + (^" A^ . ' / , ■■■ ' ■ , ■ " 

• r (N' • 'K««EG® 1 ) erT;(., 3 )=SM ■ 

C ■ C /’ f r.' 5 

£ " > K l< = 1 " ■ ' ' ■ . ' 

r/i LL tJNCCTiS (FT,rJC,N»NNf NDI F,KKfORIGINf KOUNT, SMf PHI, THETA, XK, D6Lf B’E 
: / ,WtMf »-l ,US£) 

HP sr:' = fj£. s 

C. G ' P' ■'lirPT-r AA STRIPS (IN XKI . ■ ' ■ - ■ ■' ' " 

" LL l'''4rr<fN ■ (LP,i<C,N,NN,NDIF,KK, ORIGIN, KOUNTfSMfPHI, THETA, XX, DEL, BE 
' ; ,w j.f'T AfijUSS) . 

'll I =I 
: . {*,1 ) 

J = j 

' ■ >'• " 1 = ] 

- • ' f ■ , J) = {A« { I , J) 3/DFL 

r FIRST RUN THROUGH BACKFORECASTI NG AND FORWARD RECURSION 


T Ur ' ^ 0 “ bi; 7*0 5 22 



' ■ 5 7^f fNDIFjNt 

^ ) Pkli.T '^17» (PHKI ),I = l,NP) 

’ r '-.-r •■„(■ ) PP.irsT (THM^ (1), 1 = 1, NQ) 

" f'"' ' C. ) ) PF- TfT 4-6, SM 

'Ml-, 'm,' Ce ) PF n 

“ r “ i <: . 

f +rF TG Tf^ ; 

■ , :: " = j ,1'. .T’D 

r ( ir;!), 

r. f t -rp.,. : ■ ■ ■ 

f * M ir 

tf- I'rRAjtf.r s 

I ' ' 1 T=i,K[. 

'■ I o 

^ ^ ' r J = ; f!\B 

^ . f ' , J) =• . 

M : jj=-i 
• : ' J)=AF/ (T^ J} + ?^. 
f'l JJ = !,NFi 

p r ) r : }+xy (i , jj) 

f]) ,n ) 

'* 1 t T ' )j-' . ^ 

f " G 

" F (7 3 ) Gr TO 6 £5 * ' 

r.p Tj T /G*,N'P,F:D!F,Nt, 

n. :S‘T ,1 <5,PI ,PIH,a>.,F?,»^TA,CFL , ' . 

PFTCT b^l ■ ■' , ' - . 

: . I '! 

^f'’'fTT,USS,(B Tft{ I ), J = 1,NB) ,P1[ ' 

"^^3T + .: 

F r{r"»LT*fi3 GO TO tm 

rrj!.''" r,9i ■• ■■ . .. ' : 

Of Tf tl- ■ ■' ' .■ 

c ( r ,. ■ ' 

vc-Fpr/' •^ODJFHFi (SC/=LFD AND CONSTP & INCD ) LINfARIZEO EQUATIOKS 

‘ r /. r 1 =i ,Ne 
V n c, ) )=i.G+pi 
'-MT )=l fJ)/L(I} 

r-! A<= , J = 1 , hF 

TP?;. PG, j) or TO fo: 

' r? ,j) = fA/'{i,j)/a){F3’«'0(j}3 

1 ; " n IP . ■ ■ 

!, V, r-'T F^ 

: F F.." ) C-r, TC bi.l 



' L' ' ^ ■ TKV?.^/ASfF'B,B,v fOETFR) 

^ e ) Gf- If tSH 

■ r 1 >•-.!, , 

-'TM' 

■' ''• {’ ,n=' « (’ ,1 3 

r ( “ • t 'I I . ^ , . , 

■'1 1-1,NF 

-i. C )- c.-' ' 

n. , ' = 

'■'M : ) -m 5 ( : ) + i ( F , J) ( J 3 
H( i (! 3/r (! ! 

’■ ! f TT ‘ '.r- ' , ■ ■ . 

r , < /■ ,^•F 

fT y n 

(’ 3^-fi' C( n+H(4 3 

K F r. 'i 

B LL Ui-c CKi' (BP, r.GtNjKN, NO! F,KK, ORIGIN, KDUS^TtSM, PHI, THfcTH,XXsDcLfB‘'' 

3 ^ , V' ,USS 3 

* F »USSf ) GC 1 V oil 

: V V (• C- : C" T" STING, ■ .. 

f f: 7 = i. ,N'e ' , ^ 

I F (Myr(H(I 3 )«GF.I'TA) go to 663 

"I (j. ■ ' ' „ 

f- , LL iUn.LT.-'TA , ASSUME C,ONV£RGENS* 

,:::, VC L~r 

!■'' i^r.(T )=p,f T^(j ) 

pt=p:/f2 

■: 1=1 ,NB 

a I 1 r , z " 
r>] =.PI:Vf: 

! ( G I =I tL F 
T L ( 1 ) =^B'"'r5,cn ) 

.. ■' 3=:,Kr. 

" FT 3.&Tcr'^A5 GO TO 41 

' 1 ' ■ *■ * 'iJ 

r ' /FF/P."n:T Mlf,IUM 

^ - J r „ ' 

■' : T =:1,KE 

\ • ‘ : j 

=F A (J 3* FT/«( 3=^=«'j}, ' 

t /il. U>P,NC,i%,NN,NDIF,KK, ORIGIN, KOUN 1,SM,PHI,THf*TA,XXfr/eL#B 




■ / ;. r ' , AKfUXX) 

^ ^ jUrXf (C/ TA (K) ,K=1,NB} 

~ ■; r ) --rf T. ?i 

■■ j I T ■ '.1 1, 

' i ' r 

; , - T . ,1 I' ■ 

■ r ! l.v, PJ ^ £» ) GC TO 63i 

i f I 


rr-'iM'"' P SlfiUAL VARlANCi} AMO COVARJAMCE MATRIX 

TA(I},I = l,tv’fc),PI 

•:. t ;! “ , y," 

1 K ‘ T - ■ 

' i i: ,-+;n SI DU-l. aaCD’S') 

' f ^ f UA(i ) = 

•‘r, = ; IGT!\ 

' h ’ ■ ” = ■ Q ■ ' 

tfot, (■;■!, ) l‘I;it,^^-^p+Ne+2^; 

' r:rr. 1 (.-IK + I 

j f, «AA (1)51=066, NO) 

‘'ZZ=!’- ' p> IIP 
" V ' '* =’.IS5. /FLCATCNZZ) 
f'i- 'ijSV,' R ' ; 

. •'■-< 1=1,08 . . 

r ■, .71 j=i,oe 
f i , j } = o ' 

' t - j j = i j h f., 

, ' vf , .i) = v{ I , j )+x>. (J , jj)^xx{ j, JJ) 

1'' V{ J,J) 

' rC'F* Q.I) Gf.' TC f i2 
"ill. M'' ■’■Jt'V { Vf NE,B, ,DETER) 

: F'®! ■ e ) Gf fc 673 

Of fi- I 7/ 

Tr 9 ■ . 

'■' Z Vf 1 , 3.)=r.r'./V (1,1) 

,£•/ - r I f 7 TMI ■ 

'■'! '• ] =: tC'B 

9( ’ . J) -V (1 5 J )^SVAP 

ro' " : J! = V{T.J)/SGRT(V(1,Z)’^V{J,JJ} 

) = ^:)ti-^(vi T,i ) ) 

"I- " 9 '- 5 (S (/ ) ,I = J ,NP ) ‘ ^ , 

V ;j. ' 1 " '■ '"r :n ;i\t top.-' 



M.' ■ c” To ) Sr TQ 677 

' I.' vr . c. .... ■■■ ■ . ■ 

■ C* ) to “C 6':*7 

"h , 

T' ,i- c , . 

. • , --1 ^ H' TP+Phl C J ) 

r ■ i T. , H ■•r'+ CHJ { : ) 

’H r .*.■'1' (1 ^ THi’T ) ■ , ■ ' 

■ f ^ :: «',Hric 
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li ) ■ ;■. ■' . ' .• ■ : ■ ■ 

Ai" Tt!vu? 

PRr'-r -.e^^sCGAPMA {AfU, 1 = 1, NSAPP) . 

^UGCH IIP , (GAFPA (A , I ) ,1 = 1 jNSAPP) . . ; : 

or ? i = j,r‘.SAMP ' . : 

YYYCI) =3';'-^A(4,I ) 

3 f i (i , I) 

' . (J j") =GAPP/. (/ ,1) 

/•'f': ,f.) =&APF'A{3,1) 

!1 ,•' )=GAHf' A{4, !} " ^ " 

611 ^ 

^.■i,Pn = Y , ^ 

:<:Ll P'G-PFS (HSAPP,N0PD,B3 1 \ ; 

rr 71- 1=3 ,GSAMP ., '■'-■'7:. y ^ 



V YY? • ).-i= ' " M ' ,7 ). f-' (r , , ) 

f i, T1 

. . f , ■ } ,vy ! r !. J ) 

r .- 7 1 

■ ( ' . ) .i>7-, {/,,.) . . 

,,]•■ '■ ■' . 

^ • if' •',! C--'^ , , 

' ’ t i. ^ ‘ ‘ 7 ^ 7' 1 P f 1 7 p £ ) 

•i, 

rr r!*.'- Y- J h BU) *G AMF./. ( 1 , 1 ) « B ( P }=^GA 1 f Z f 

' ’ J- ( : J ) 

I- n ii': 

'••’F-'.i'’" ■■ u' i = 

fffTB nfi=]pr r-iSAF) ■ ■ 

Af l.:f ■ ■ . » ' 7 :■■■ •:: : / 7: 

YYYf fppj) ^ ■ ■ ' . ■ ■ 

Y .'MT ,1 ) =r,AMH^ dpi) 

V'v ( J 5 t£p T ) 

;•’/ (I ) *GAMf^A(?- y I ) 
v-'f I ,/ I =G *MWA (4t I) 

V' (I t*' I »CA^MA(5s I) 

6U CCAITifUt' ■ . : 

WC.K,D*5 ■ ■ . 

C^LL R*' GRf,S{NSAWP,NCPD,B) 

DC 617 I=i,NSAMP 

YYY( I )=eFTaj 6,n«GAf'RA<6»n 

vy{I #1 ) 

?'> {I ,?.}»GAMHM2,1) 

r.yAi 

vyif p/. ) (A , n 

YXCI ,7) »GA^'^^A(5 ,1) 

617 COKTIWC ' , ■ ■ 

MOROaf 

C'ALL RFGRtS { KSA?'Pf ftCRC, B ) 

PG W I«l,NSAMP 

GAKiHA(-.> fl )ae'.TA{Y^p I ) " B ( 1 } *GA.wMA( If I)'»B{ 2 l*GAH>^A( 2, I J-BC 3 3» 

?!'J*"8W )*GARiYA(4f I)-bC5J>i'GAMMAf5,I) 

616 CCNTIKMJF 

PRIWT (GAI^MACGf i ) 1 1 =l ,NS AMP ) 

P'UNCH 77s , (GAMHa (6, 1 ) , I = 1,MSAMP) 

' do e2Y I = lf^,SAMP • ' . ■ 

tYV{ T 5*bf=’TA Ccfl) '' . 

' ^ ,'^Xf I fl, I'^GAMKAdpI) i '' 

' ' XXd »2) =GAFMA(?5 n 

: I'tr: j3)=om-n3ti) 

1 tXf ; 7*7) aCr^jyf'/ (6 , 1 > 

'.'7;' ' tni tA ) aG7-*1Ri.d6f i> 

/.62f'-CC‘7;-: d!'.-?-: 

: rfp. >n:kt'-\ - ' • 

» [ 1 ': 


r ^ ^ rr ij[ s i kOf'Df B ) 

Blip 


if ^ . s.’t-r*#. 7l£w 5^* ‘ 3|e». ;:^.j». :^w 5 %-b s ^51 ^ - 3^;.. 4= - 

5ljf: j MiTJfjr grC'Pt S ■■ ■ •' ■ 

“h;^ ':UfFr:UTTr.; JS C^LL^D EY THE MAIN PROGRAM FOR R FGRFrSSIlif: AK’tLYS ISo 
’'■/r^ i,1!Pe»'R OF IH’c DATA POINTS. 

(i-:i ’H. FO. UF INDeP3NDfiHT VARIABLtS- 

5j, n *1* - Z'* , y' ‘ ^ . *11 ^.tur f it :4C«* J^ir 3^«3 ^r-j ^wp^S^as, 


^ Uf.i- f. LV 1 nr R*:-&RLS {KSAPP,NORC,BeTA) 

0 ) '^: pH Y < 250 ) t Y C F 5r.. , I U ) t KTR P( 1 'i , 3 53 ) f YEST ( 35 T ) t KTR PK { li*F, 1 ~ ) 

i.JM'' ;■ MCANd* }fXTRPY{lO},BETA{U' J,8{lC)t6B(l.l, J),SUMX(i*)), 

’Ml) 
r i..F.v'>rF-‘ YjX 
;■: /L Y’ SOR 

1 Ts-'f ■ , ' ■ , 

Of. L^^ljfsCRD 

: (1) 

on ?,r 1=1 tfJSAMP : .. , 

sup:' (L r»S!.iFX(Li+i;(i,L) 

rr.MlMJF 

'’MFA!SfL)*SUMY{L)/FLCAl (NSAMP) 

2/. Cf'NYjMCjr 

SUFY«<d5‘ . . ; , : ■■■■ ■: ■ 'o;: ';0: ' 

no 4fi I =1 fNSAMP , ' . 

tli|«Y='‘UfiY+Y fl 5 ■ ' 

4u COKTIMUL 

YMf AN*SUMY/FLOAT (FSAMP ) 

PRINT lC»5,(XMe-AN(L), L=l,NORD>,YMEAN 

Dt. 320 l*I,NORD '■' ' ■: 

DO 3Do I=1,NSAMP 
:':{|»L) = K(I,IHXMFAN(L) 

3?><'< COhUMH ■■ . V -n: : vO;:: ; 0' -: 

DO 5r l=l,NSAMP 
YU? »Y( D-YMfAN 
C< CCiKTINUn 

HP r.0F'YTf-'uF ■ 

OO'i'tift J=l, yNSAMP 

DO t€ L-l. fNORD - /:'■ 

xTRpaf ij=xOfL) ■ ^ -r;:- 

&f: cr.'N^JAUE ' , •■ ■ 

00' •’C I=J. rNCRD ■ ^ 

.0G'7f? L=rVNCRC ■ ' ' . i 

KTRPXn sLl'=0»0 ^ 

.'POnTFi K=1,NSAFP -■ 

'KTRPXO fL)=XTRPX(I,t)+XTRPri,KI*KCK».i''i ,■ 



■ i . 'f'.ijr 

" r f Q* 1 ) Gt Tf ir 

■' ' 

YY=- c 

T’f I =1 siiS^PP 

"■ YV-SYY+Y( I ) ^^2 ■ 

r-f Gr L = l,^s(:RD 
^fDr’^TPP,. (L,L) 

CGaTIMJ;: 

l''5tSVY» {S(L),L = 1,N0RD) 

fi T=l,GSAf«P 

Y(?)=YC a}/SCRT {SYY) 
i‘ c I f’"'' Tf 'I'r 

n( 1. ?! L = l,|.Ci'.D 

PU 5 2'' 

■f JjL):^; {J,L)/SGR''-{c.(L) ) 

12. 

Gf- YC 

13 PYiL '^PTJf^iv «XTRPY,PiCRDfB8fjT,DErRM) 
r-l; lyr' L=I sNCRD 
,'^'TpPYfL)=t .P 
MP r* r Tjrufr 

DC 5. ‘ji. I = :tJ^CRD 

Dr'lSC, J=l,RiSAMP 

''TRPYfl 3*KTRPY(I 34XTRP(ItJ)=«'Y{ J) 

15t rt-KTINUE 

Df' IBt:* L*|»f^ORO 
BCDaCcP 

K*r rfiPX7f-yp 

DC! 17./ L=ltNDRO 

Drja7lft".-‘J«lfNCRD ,■ 

B a)»e(L> + XTRPXCL, JJ’^XTRPYt J) 
pf rriMiR'uE 
PRINT 106 

PMffT 1«,5,{E(1 ),I=l,NOBD) 

Dft J PG L=liKC'RO 
BETA CL) =Ba}*SQRTfSYY/S(L } 3 
if'O oC'O'^J.r'ue- ■■ ■ 

■■ PRINT 11-7 ■ ■ . 

•PRINT Itv^ffBFTAUJyL^lfNCRO} • 'P- 
suNYY!=r}*n' 

DC- -251)’ l^lf^SAPP 
SUFYY^SIJYYY+YU ) -■ - 

?l<- CdMiG'Uf • - 

■ > " YKNDR*eUMYY/FLOATUJ3AHP3 P ■ 

OPR JNt 2"9ifYRN0R 
■Or SSRrc»‘J*CrP 
1, DO l=»l 



'■ ^ ■; " J = ! ir-.CRD 

J7(i )+B{ If J) 

; " < i " ij-" 

:?-'i-G+ 1 YFSTd )' Yf^NCR)=«‘*2 

i"' ; T'iini:” 

■ f' :■•:• : =’ ,k?j-fp 

r r r' ■;!■=:: , 4 jr,+ {Y (I )■ YMNrR) 5 )= 5 it 2 

' ; ’ll ri> .'lf-^:i^^•'SSR'■G 
’■! F*" !■■■=.' ‘S/'MP- 1 
"Drr f-G^rTi}- 0 
’T d I’ F>=!'inFlCT-NDFFr& 

':QP = :.GB’'-G/FL0AT{MCFRFG3 

YAP" r. ^Rr S!UL/PLr^T(NDFRFS ) 

FVAUJ' =R;»'S0P/VAf'. JAN 
f' I- : 1 }. ! F 
PPJl ■’■ 1 G 

Pf I'f'.r’ 1 ? f SSK ANfNDFTOI 
"T'Tf,? If 2,SSftFG,NDFREG,MeSQR,FVAlUE 
'•■l-l'd ’F^fPPRJDUfNDFRFSfV/.RIAN 
l-;d FI P^r./>*d2X,*T( TALICCRRFCTFO FOR MEAN) *, 2Xt £16.7, 2X, I 5//) 

1 ''A FGf-.M/.''{l.?X,*RfGR’"5Sir!N‘!',5X,*' 16. 8, 13*( , 15, 5X, S 16. 8, 5X, t: 16. 8,5X,cl6.e 
1//) 

ir ' F OR M A “ f 1 » X , PF S I DUA L=i-',!--X,£16.8,lUK,I5,5X,El6.8) 

i>' FfPM4‘’'( i5X,=!‘SrjUFCP=^,8X,=>=SUM CF SQUA^FS’i'tSX, ^DEGREES OF FREIOOH#,5X 
S€UAPF^,S Xt^F VALUF=<'//) 

Iff Fr-P.MAT(1X,1< L13.5) 

IC 6 ■ FCPMATn HF OLD PARaMCTFRS are*) ' '' 

r;7 ACTUAL PARAMETFPS OF THF MODEL ARE*) 

lf« FORMAT t5i,^y,>KANAL YS IS OF VARIANCE TABLE=<') 

If 9 FrPMATdXfYtHO X TRANSPOSE X MATRIX’S') 

25r F(,!RMAT(lX,*Mf-/N CF NORMALIZED YCI) SERIES =’^=,£16.8) 
r'-TURF 
TOC 

it*' «> :^rai 3$: •> 4i u 4; IR ’S' w ’^*>'>’1'*. 't'. 


C SUBrC'*dTN£ CXRrLA • •' ' ■ ■?■- '■ '' 

this’ SIBRDUTINE CALCULATES CRCSS-CORRE tATlCN FUNCTION BETWEEN ALL PAIRS'^ 
THF V'TTf.'RSo 

w #,*s, :|c*w;tC(w 

MlAG ... .KAXIPUM NC. OF LAGS* : ■ 

Rfd) «,oo CLRK'LATTCN COEFFICIENT. , l.Lf'G'/ ^tv, 

, WPfl) CCVARIAI-XF COEFFICIENT. 

,, ■ Sl’PP CUT INF. CORFLA { ASAMR, , ''T '• 

D I Mr K ST CN C R j ,WP I «P ) •• Tf , ■';: ' '* 

DJ MF t\$I tm . X f 12,353 ) , SP( FP C I5& ) 


!-;f = I,KVCTP 

" Y = j. jf.VCTR 
i' f] ) ^ - 

ni 

TMi )= . - 

3 ^ T iHZAhi- 

" ? C 0CI)+XC!K,I ) 

(1 }+> ( IY,1 )»*2 

Ffi: )-Fpn )+>■(! Y, I ) 

1 r : , !’t:!:i!; 

Y1rf^LAG + X 
ri 7 e 

J = J'”7 

Krf.'S .» f^P-T +J 
•■MI) =^T?( J)' X{ 

?P (I )=.FrM!- 1 j.-x< IXf J)’^*2 
Pp !I P{ J)->:tIY,K) 

■ "vpa )=GPCr-1 3»X{!Y,K)**2 

2* r.oTir-'UF ■ ■ 

Ml/'G1=HUC+1 

Dl ;M f=l,f1LAGl . .. 

';:NINP=r:SAFP« l+l ■'■■'■■':- 

CFO . . . - ■ ■ 

PC 26 J=lf^KINP ■ ■ ■ ■. ■ ■. . 

K2 * J + J" 1 ' ■ ■■■■■' 

26 CPCU»CP(I)C{X{IX,K2)*X(IY,J)) 

WP(I ?«CP{|)/FLOAT(NKINP) 

M;(X’;®'=LCST (6 MINP)*CP{ n-FPin^lPd ) 

RLCNl-SQP'^t {FLCAT(NHJNP)=4'GP( I ) )“FP{I )**2) 

RDB)2*SQR1C jFLCAT(6‘HIMP)=S'SPn ) HTPd )**2) ' 

%(' RP{! )s«RNUF/ (PD6Fa=^Rr<EN2) 

I, MIX-IY) 2r’ 2, 262, 2. 2 
2f>2 PPINT ]£r,IX,lY 

^ . GMTC 2i 4 ‘ 

20? PRINT <3'5,!X,IY . 

2r 4 C6f*‘TTI'iU2 

99' FOPK/iTUXjY^LiTO CORkELaTICN FUNCTION BETWEEN Xs*»I4,*4N0 Y=*,I4} 
Jf'i FOPHAT ClYj^CROSS PrppFtATION FUNCTION BETWEEN K=*,I4,*AND Y=4,I4) 
PRIN* 6&6,(PP(I),I=1,MLAG) 

660 rtPHATf 1^,1 6F8o4 ) 

CH|SQT~f'o‘. ' ^ ■ . ■ 

CHTSC2=i,)«< 

' t'HlSe' =C'.f . . ■ 

OO fir fC.Hl = !, ,?,f J' . ^ _ 

lYsTfHJ+J , . ^ . 

'i’F'' (|%L£'.13) CH!SQl = CHli«i.'FRP'llT|''I«»lY}**2 
'If ClJolFoPS) CHISQ2»CHISC:2r+RPftl.tX»iyi4*2 
't'F:MnoL'“«31 ) CHlSQ3»CHlSQS+'.t^PtI.I»W'flYJ4’i‘2 



: ' •: "iMi- ■ : ■ ■ ' , ' . ' . ■ ■■ 

f-iU(J3=^FL0/iT 5 tsSAFP) : 

Ui’ : :g2=«'FL! ■/*"■{ •'^SAFP) 

•'I'':: ii.I so:- «FLt AT( f.SAFPI 

rrMi-jsr ?r^e) 

H'. T e6,UHJ:,TCHli,TCHIi ' ■ 

f }r,=5-<‘.H-’SCUAF^' statistics FR3M lag 1 T3 l2=>«',E2t’. 
statistics from lag i to 2a==f, / ix, * chi»souafh ST*”ir 

lag 1 TO :.r = =^*E2ii4.8) 

( "T-'HI. ■ . : : ■ 

‘ , ■' ■ , ■ 

rr. ■ ' ' 

* 7/ « r S^' 15 ^ ^ ^ « ^iK*- 3^ w !l^«» 5^i» 3fe»* :^m sft m s") 'V * '1 \s> « 



